Time Dependent Hadronic Modeling of FSRQs

C. Diltz11affiliation: Astrophysical Institute, Department of Physics and Astronomy,
Ohio University, Athens, OH 45701, USA
, & M. Böttcher22affiliation: Centre for Space Research, North-West University, Potchefstroom, 2520, South Africa 11affiliation: Astrophysical Institute, Department of Physics and Astronomy,
Ohio University, Athens, OH 45701, USA
, G. Fossati33affiliation: Department of Physics and Astronomy, Rice University, Houston, TX 77251, USA
###### Abstract

We introduce a new time-dependent lepto-hadronic model for blazar emission that takes into account the radiation emitted by secondary particles, such as pions and muons, from photo hadronic interactions. Starting from a baseline parameter set guided by a fit to the spectral energy distribution of the blazar 3C 279, we perform a parameter study to investigate the effects of perturbations of the input parameters to mimic different flaring events to study the resulting lightcurves in the optical, X-ray, high energy (HE:  MeV) and very-high-energy (VHE:  GeV) -rays as well as the neutrino emission associated with charged-pion and muon decay. We find that flaring events from an increase in the efficiency of Fermi II acceleration will produce a positive correlation between all bandpasses and a marked plateau in the HE -ray lightcurve. We also predict a distinctive dip in the HE lightcurve for perturbations caused by a change in the proton injection spectral index. These plateaus / dips could be a tell tale signature of hadronic models for perturbations that lead to more efficient acceleration of high energy protons in parameter regimes where pion and muon synchrotron emission is non-negligible.

galaxies: active — galaxies: jets — gamma-rays: galaxies — radiation mechanisms: non-thermal — relativistic processes
slugcomment: Accepted for Publication in The Astrophysical Journal

## 1 Introduction

Blazars are a subcategory of active galactic nuclei that can be divided into two general classes, namely BL Lac objects and Flat Spectrum Radio Quasars (FSRQs). They are typically characterized by their rapid variability, superluminal motion and their extreme luminosities, often dominated by their -ray emission. These features are considered to be the result of beamed emission from a relativistic jet oriented at a small angle with respect to the line of sight (Urry, 1998). The broadband spectral energy distributions (SEDs) of blazars can be typically characterized by two broadband, nonthermal peaks that span from the radio to UV wavelengths and from X-rays to high energy -rays. It is generally accepted that the first spectral component from radio to UV wavelengths is the result of synchrotron radiation of electrons/positrons in a magnetic field. For the origin of the second broadband peak, two different paradigms are often invoked, collectively referred to as leptonic and hadronic models (Böttcher et al., 2012). In the leptonic scenario, the high-energy (X-ray – -ray) emission is due to inverse Compton scattering of low-energy photons off the same electrons/positrons. The low-energy target photon fields can be the synchrotron photons within the emission region (SSC = synchrotron self Compton), or external photons (EC = external Compton), which can include the accretion disk (Dermer et al., 1992; Dermer & Schlickeiser, 1993), the broad line region (BLR; Sikora et al., 1994; Blandford & Levinson, 1995), or infra-red emitting, warm dust (Blazejowski et al., 2000). Leptonic models have been quite successful in explaining many features in the SEDs and light curves of blazars.

Hadronic models have also had success in modeling of blazar emission (e.g., Mannheim & Biermann, 1992; Mannheim, 1993; Mastichiadis & Kirk, 1995, 2005). The hadronic model suggests that the high-energy emission originates from the synchrotron emission from a ultrarelativistic protons. The relativistic protons interact with the radiation fields within the emission region, producing high energy pions, which then decay to produce muons, electrons, positrons, and neutrinos. The pions and their decay products emit their own radiation (primarily synchrotron) which adds to the broadband spectral components in the SEDs of blazars (Aharonian et al., 2000; Mücke et al., 2003).

As shown in Böttcher et al. (2013), leptonic and hadronic models are generally successful in reproducing the SEDs of many -ray blazars. Therefore, one needs additional diagnostics to distinguish which type of model is most applicable to blazars. The most obvious difference is the production of TeV – PeV neutrinos produced only in hadronic models. Additionally, due to the vastly different acceleration and cooling time scales expected for electrons/positrons vs. protons, one also expects substantially different variability patterns predicted by the two types of models. This latter aspect is being studied in detail in this paper. Note that an alternative discriminant may be the polarization of the high-energy (X-ray – -ray) emission of blazars, as discussed in Zhang & Böttcher (2013).

Determining the underlying shapes of the particle distributions that give rise to the broadband spectral components, is critical to understanding the physics of particle acceleration and cooling in AGN jets. Simple power law and broken power law proton distributions can be produced through diffusive shock acceleration when incorporating radiative losses, and such distributions have been invoked in hadronic models to explain the emission and subsequent particle cascades that produce high energy -rays in blazars (e.g., Mücke & Protheroe, 2001; Böttcher et al., 2013). Second order Fermi acceleration is a viable mechanism for producing log-parabola-shaped, curved spectra. The curvature of the spectra can give clues to the parameters governing the Fermi II acceleration mechanisms (e.g., Schlickeiser, 1984a, 2002).

Recently, a time dependent hadronic model that considered a Fokker-Planck equation with the incorporation of radiative losses, second order Fermi acceleration and the emission produced from the final decay products of the photo-hadronic interactions was utilized to explain blazar emission (Weidinger & Spanier, 2015). The production rates of final decay products were derived by analytical parametrizations of the energy distributions for the neutrino, electron/positron and photon distributions from the interactions of relativistic protons with the photon fields, Kelner & Aharonian (2008). However, in order for this approach to be viable, the synchrotron cooling time scales of the intermediate decay products (pions and muons) must be significantly longer than their decay time scale (in the co-moving frame of the emission region), which restricts the combination of maximum proton Lorentz factor, , and magnetic field to  G (Böttcher et al., 2013). If blazar jets are the sites of the acceleration of ultra-high-energy cosmic rays ( eV), then such models are only applicable in the range of magnetic fields of  G, substantially lower than usually found in hadronic modeling of blazars. For higher magnetic fields or maximum proton energies, the synchrotron emission from muons and possibly also charged pions becomes non-negligible. At the time of writing, no time-dependent hadronic model has been published which incorporates the radiation emitted by the pions and muons produced in photo-hadronic interactions.

In this paper, we describe a new, time dependent hadronic model that considers the radiation emitted by all secondary products and incorporates Fermi acceleration and self-consistent radiative losses for all particle species (including photo-pion production losses for protons). We use this new code to provide a fit to the average SED of the FSRQ 3C 279 to determine a baseline parameter set. We then apply a Gaussian perturbation to several input parameters (specifically, the magnetic field, the proton injection luminosity, the Fermi-II acceleration time scale, and the proton injection spectral index) in order to study the resulting light curves in the optical, X-rays, HE and VHE -rays as well as neutrinos in the energy range of sensitivity of IceCube, and study the characteristic variability patterns and cross-correlations between the lightcurves in the various bandpasses. We describe the model setup in §2; we present a model fit to the SED of the blazar 3C 279, using an equilibrium solution of our code, in §3; we then study the light curves resulting from the Gaussian perturbation of the selected parameters in §4, and compute the cross-correlation functions between the various light curves, in §5; we present a summary and brief discussion in §6.

## 2 Model Geometry and the Time Evolution of the Particle Spectra:

### 2.1 General Considerations:

In this section, we describe the assumptions made in our model, the features of the Fokker-Planck equations used for the evolution of each particle distribution, and the radiative components that they produce. We assume a homogenous, one zone model, where a population of ultra-relativistic protons is continuously injected into a spherical region of size , moving along the jet with a bulk Lorentz factor , embedded in a homogeneous, randomly oriented magnetic field of strength B. The size of the emission region is set by the observed minimum variability time scale, , through

 R=c⋅Δtvar⋅δ1+z (1)

where is the redshift to the source and is the Doppler factor. The time evolution of the proton distribution is described by a okker-Planck equation that incorporates cooling due to synchrotron radiation and pion production. The proton distribution interacts with the photon fields and generates relativistic pions. The pions subsequently decay to produce muons and muon neutrinos. The muons themselves also decay to produce electrons, positrons and muon and electron neutrinos. The evolution of each of the particle populations is described by their own Fokker-Planck equation.

We assume that the proton injection spectrum takes the form of a power law distribution:

 Qp(γ)=Qp,0γ−qpH(γ;γp,min,γp,max) (2)

where denotes the Heaviside function defined by if and otherwise. The normalization factor for the proton injection spectrum is constrained through the proton injection luminosity, :

 Qp,0=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩Lp,injVbmpc22−qγ2−qp,max−γ2−qp,minif q≠2,Lp,injVbmpc2ln(γp,maxγp,min)if q=2. (3)

where denotes the comoving volume of the emission region and denotes the rest mass of the proton.

All particles can be accelerated or decelerated by gyroresonant interactions with magnetohydrodynamic waves. This interaction causes the particle distribution to diffuse in energy, pushing particles to higher and lower energies. If the energy density of the plasma waves is small compared to the energy density of the magnetic field (quasi-linear approximation), then the diffusion coefficient governing the momentum diffusion mentioned above, takes the form of a power-law,

 D(γ)=K⋅γp (4)

where the proportionality constant is set by the shock velocity, , and the Alfvén velocity, :

 K=1(a+2)tacc (5)

where . We invoke a diffusion coefficient with a spectral index of (hard sphere scattering). This allows the acceleration time scale to be independent of energy.

### 2.2 Pion Production Templates:

The protons also interact with the photon fields and produce neutral and charged pions. The total proton-photon cross section is divided into separate components, corresponding to different reaction channels through which the pions are produced: direct resonances (such as the resonance), higher resonances, direct single-pion production and multi-pion production. We use the prescription developed by Hümmer et al. (2010) for the photo production rate of pions:

 Qb(Eb)=∫∞EbdEpEpNp(Ep)∫∞ϵthmpc22Epdϵnγ(ϵ)Rb(x,y) (6)

where denotes the proton distribution, denotes the photon field that the protons interact with as a function of the normalized photon energy , and (corresponding to an energy of 150 MeV) represents the threshold below which the cross sections are zero. The dimensionless variables and are given by

 x=EbEp (7)
 y=Epϵmpc2 (8)

The response function in the photo production rate of pions is given by

 Rb(x,y)=∑ITRIT(x,y)=∑IT12y2∫2yϵthdϵrϵrσIT(ϵr)MITb(ϵr)δ(x−χIT(ϵr)) (9)

and is summed over all interaction channels that make up the proton-photon cross section as a function of photon energy in the parent nucleus rest frame, . The functions and represent the multiplicity of daughter particles and the mean energy fraction that is deposited into the daughter particles for a given interaction channel, respectively. Evaluating these integrals turns out to be very cumbersome. Therefore, Hümmer et al. (2010) suggest a simplified prescription in which the interactions are split up into separate components that take into account the resonances, direct production and multi-production channels and which assumes that the multiplicity and deposited mean energy fractions are independent of the interaction energy, . The response function then simplifies to:

 RIT(x,y)=δ(x−χIT)MITbfIT(y) (10)

with

 fIT(y)=12y2∫2yϵthdϵrϵrσIT(ϵr) (11)

With the simplified response function, the photo-production rate of pions can then be written in the more compact form:

 QITb=Np(EbχIT)⋅mpc2Eb∫∞ϵth/2dynγ(mpc2yχITEb)MITbfIT(y) (12)

This single integral is easy to evaluate numerically. The photo-production rate of pions now depends on the response function, , the multiplicities, , and the mean energy fraction deposited into the secondary particles, , for the dominant interaction types. The values of the multiplicities and the mean energy deposited for the resonance, direct production and multi-pion production as well as the response functions used, are tabulated in Hümmer et al. (2010). With the formalism adopted for the pion production rates, the cooling time scale for the proton distribution from the production of pions is given by:

 t−1cool(Ep)=∑ITMITpΓIT(Ep)KIT (13)

where is the inelasticity and the interaction rate of a given interaction type, given by:

 ΓIT(Ep)=∫∞ϵthmpc22Epdϵnγ(ϵ)fIT(Epϵmpc2) (14)

Using this formalism, the rate at which primary protons are lost due to conversion into neutrons can be given by the expression:

 t−1esc,n(Ep)=∑IT,p′≠pMITp′ΓIT(Ep) (15)

where denotes the new nucleon created in the photohadronic reaction. This formalism also allows the energy-loss term in the proton Fokker-Planck equation due to pion production to be defined as:

 ˙γpγ=−Ep⋅t−1cool(Ep)=−Ep⋅∑ITMITpΓIT(Ep)KIT (16)

Given that the radiative cooling time scales for protons can be longer than the typical dynamical time scale of the expansion of the emission region, we include adiabatic losses in our model. Assuming a conical jet with an opening angle of , the adiabatic cooling rate is . The full proton Fokker-Planck equation that incorporates radiative losses due to synchrotron, adiabatic processes, pion production, neutron production as well as stochastic diffusion through the interaction MHD waves reads (Schlickeiser, 2002):

 ∂np(γ,t)∂t=∂∂γ[γ2(a+2)tacc∂np(γ,t)∂γ]−∂∂γ(˙γ⋅np(γ,t))+Qp(γ,t)−np(γ,t)tesc−np(γ,t)tesc,n (17)

where denotes the dynamical escape time scale for the protons which we parameterize as a multiple of the light crossing time: where . The value of is kept as a free parameter. The term denotes the combined loss rates on the proton distribution due to adiabatic, synchrotron and pion production processes. The synchrotron loss rate is given by

where is the energy density of the magnetic field.

### 2.3 Pion and Muon Evolution

The decay time scale of neutral pions (in the pion rest frame) is only  s, and they are not subject to synchrotron losses. Therefore, in our code, neutral pions are assumed to decay instantaneously and so, no Fokker-Planck equation needs to be solved. The charged pions, however, have a significantly longer half life ( s in the pion rest frame), so a separate Fokker-Planck equation has to be considered for the charged pions produced in proton-photon interactions:

The only major difference comes from the loss term due to the decay of charged pions with a characteristic timescale . If the Lorentz factors of the pions are large enough, the decay timescale could be of the order of or even larger than the pion synchrotron cooling time scale.

Charged pions decay to produce muons through the following channels:

 π+→μ++νμ (20) π−→μ−+¯νμ (21)

The decay term in the pion Fokker-Planck equation (last term in Eq. 19) serves as the injection function for the muon Fokker-Planck equation. The muons then follow their own Fokker-Planck equation that incorporates loss terms due to synchrotron processes and diffusive acceleration. The muons can then decay through the following channels:

 μ+→e++νe+¯νμ (22) μ−→e−+¯νe+νμ (23)

to produce separate distributions of electrons, positrons and electron and muon neutrinos. In total, we have Fokker-Planck equations for the proton, electron/positron, muon, pion, and neutrino distributions that are all coupled to each other and represent all particle populations within the emission region. Note that the electron/positron Fokker-Planck equation contains an additional injection term due to absorption and pair production, which allows us to properly follow the evolution of ultra-high-energy gamma-ray induced pair cascades (see, e.g., Böttcher et al., 2013).

Once we know the individual particle spectra, we can compute their radiative output (primarily due to synchrotron emission) at any given time step. The synchrotron emission coefficient for a distribution of charged particles within a tangled magnetic field is given by:

 ji,syn(ν,t)=14π∫∞0dγni(γ,t)⋅Pi,ν(γ) (24)

where the term denotes the synchrotron power per unit frequency produced by a single charged particle of species , and can be well approximated by (Böttcher et al., 2012):

 Pi,ν(γ)=32πc9Γ(4/3)r2e(memi)2⋅uBγ2ν1/3ν4/3ce−ν/νc (25)

where denotes the energy density of the magnetic field, the classical electron radius and the mass of a particle of species . The critical synchrotron frequency , is given by

 νc=4.2×106⋅B(G)⋅(me/mi)⋅γ2 Hz (26)

The synchrotron spectrum represents one component of the combined photon field that also includes the synchrotron-self-Compton and external-Compton radiation of the electrons/positrons (Compton emission from the heavier particle species is strongly suppressed due to their much higher masses) as well as the radiation fields produced by the decay of neutral pions. We then solve a separate evolution equation for the combined photon field (Diltz & Böttcher, 2014):

 ∂nph(ν,t)∂t=4πhν⋅∑kjk,ν(t)−nph(ν,t)⋅(1tesc,ph+1tabs) (27)

where the sum is over the all radiation mechanisms and all particle species, denotes the photon escape time scale and denotes the absorption time scale due to synchrotron-self-absorption by electrons and absorption. The absorption time scale can be defined through the opacity as

 tabs=Rc⋅(τSSA+τγγ) (28)

where and denote the synchrotron-self-absorption and absorption opacities.

We utilize the head-on approximation to simplify the differential scattering Compton cross section. Using the head-on approximation and assuming that the electron distribution and synchrotron photon fields are isotropic in the comoving frame of the emission region, the comoving SSC emission coefficient is calculated following Jones (1968). The incorporation of the external radiation fields is implemented the same way as in our previous paper on a time-dependent leptonic model (Diltz & Böttcher, 2014). We compute the absorption opacity using the prescription of Dermer & Menon (2009), and the pair production spectrum as given in Böttcher & Schlickeiser (1997). The produced pair spectrum is added to the solution of the electron/positron Fokker-Planck equation at every time step.

With the combined photon field at every time step, the components that represent the broadband spectral energy distribution are then found through:

 νFobsν(νobs,tobs)=h⋅ν2⋅nph(ν,t)⋅δ4⋅Vb4πd2L⋅tesc,ph (29)

with and .

### 2.5 Neutrino Emission:

Our code also takes into account the production rates of electron and muon neutrinos generated in muon and pion decays following photo hadronic interactions. The neutrino production rate depends on the number of charged pions that decay within a given time, , which is given by the decay term in the pion Fokker-Planck equation:

 Dπ(Eπ)=nπ(γ,t)t′decayγ (30)

With this pion decay rate, the neutrino production rate can be calculated as

 Qν(Eν)=∫∞Eν(1−rM)−1dEπEπ⋅Dπ(Eπ)1−rM (31)

where .

The rate of muon decays is governed by the muon Fokker-Planck equation. The calculation of the spectrum of neutrinos generated by the decay of charged muons is more difficult than in the case of pion decay, since the system is a three body decay. We follow the procedure derived in Barr et al. (1988) to find the neutrino production rate for the three-body decay of muons:

 Qν(Eν)=∫∞EνdEμDμ(Eμ)dndEν (32)

Using the dimensionless scalar variable , we can recast equation 32 into the form:

 Qν(Eν)=∫10dmDμ(Eν/m)m⋅dndm (33)

where represents the neutrino production rate in the laboratory frame in terms of the dimensionless variable . Assuming that the neutrinos are traveling relativistically, we can cast the neutrino distribution function into the following form (Barr et al., 1988):

 dndm=g0(m)+g1(m) (34)

The scalar functions and are listed in Table 1 and describe the laboratory-frame distributions of the neutrinos in the relativistic limit. Once we have computed the neutrino production rates within the emission region, we determine the expected fluxes here on Earth and integrate over the IceCube sensitity range in order to determine the expected number of detectable neutrinos.

## 3 Application to the FSRQ 3C279

3C 279 was the first -ray blazar discovered using the Compton Gamma Ray Observatory and has been the target of several multifrequency campaigns (e.g., Maraschi et al., 1994; Hartman et al., 1996, 2001; Ballo et al., 2002). 3C 279 has been classified as a flat spectrum radio quasar given the location of its synchrotron peak in the infrared. Many observational properties of 3C 279 have been well measured, including the accretion disk luminosity (Hartman et al., 2001), the bolometric luminosity of the broad line region (Xie et al., 2008), the minimum variability time scale (Böttcher et al., 2007) and the apparent superluminal speed of relativistic jet components (Hovatta et al., 2009). 3C 279 is one of only three FSRQs detected in VHE -rays by ground-based Cherenkov Telescope facilities. Specifically, 3C 279 was detected by the Major Atmospheric Gamma-Ray Imaging Cherenkov (MAGIC) Telescope during an exceptional -ray flaring state in 2006 (Albert et al., 2008). Böttcher et al. (2009) have pointed out that this VHE detection, in combination with the rest of the SED and known variability properties of 3C 279, presents a severe challenge to single-zone leptonic jet models, and suggest a hadronic scenario as a viable alternative. For this reason, We choose this well-known blazar as a representative of -ray bright blazars in which hadronic processes might be important.

We performed a parameter study to provide a rough fit to the average SED of 3C279 (as presented in Abdo et al., 2010), by running our time-dependent leptohadronic model code with time-independent input parameters and waiting for all particle and photon spectrum solutions to relax to an equilibrium. To obtain these equilibrium solutions, we set the size of the time step in our code initially to  s. This time step size is considerably longer than the time scales for all loss terms, acceleration terms and escape terms in all particle Fokker-Planck equations. The implicit Crank-Nichelson scheme used to numerically solve the Fokker-Planck equations guarantees that the simulation converges to a stable solution after a few time steps.

Given the number of input parameters in our model (see Table 2), it is important to independently constrain as many parameters as we can from observational data. For 3C279, we have the following observational parameters (see Böttcher et al., 2013, for references to the observational data): , (the apparent transverse velocity of individual jet components, normalized to the speed of light),  d,  erg s, and  erg s. The superluminal motion speed sets a lower limit to the bulk Lorentz factor, . The observing angle is set by using the relation so that . From the variability time scale, we can constrain the location of the emission region along the jet, . With the luminosity of the broad line region, we can determine the characteristic size of the BLR using the luminosity-radius relation (Bentz et al., 2013). The mass of the supermassive black hole in 3C 279 is constrained through the measured bolometric luminosity of the broad line region and is found to be (Woo et al., 2002). With the mass of the black hole and the accretion disk luminosity, we can then constrain the Eddington ratio for the accretion disk emission. We approximate the BLR spectrum as an isotropic (in the AGN rest frame) thermal blackbody with a characteristic temperature of , which has been shown by Böttcher et al. (2013) to yield Compton spectra that are virtually indistinguishable from spectra using more complicated BLR reprocessing geometries.

Within the parameter constraints listed above, we perform a ”fit by eye” to find suitable values for the remaining parameters. In the context of most hadronic modeling, the X-ray to soft and intermediate -ray emission from FSRQs can be best explained by proton synchrotron radiation. Thus, the X-ray through HE -ray spectrum informs our choice of the proton injection luminosity, spectral index, and maximum proton energy. The VHE -ray emission detected by MAGIC (Albert et al., 2008) appears to constitute a separate radiation component beyond the Fermi-LAT -ray spectrum, and we here suggest that this component may be provided by muon- and pion-synchrotron radiation, which our code is uniquely able to handle in a time-dependent fashion. By chosing , our simulations will be in a parameter regime in which muon and pion synchrotron is expected to make a significant contribution to the -ray emission. A full list of parameters which yield a satisfactory representation of the SED of 3C 279, is given in Table 2.

With this set of baseline parameters, the broadband SED of 3C 279 can be reproduced quite well, as shown in Figure 1. The infrared to UV portion of the SED is fitted by synchrotron radiation from primary electrons. The X-ray to GeV -ray emission in our model SED is dominated by proton synchrotron radiation. The VHE -ray spectrum, as measured by MAGIC, is best explained by a combination of synchrotron radiation from the primary protons and secondary muons and pions generated via photo-pion production. We note that the proton synchrotron component slightly overshoots the Fermi data points. This is reasonable since the Fermi-LAT spectrum represents a long-term averaged high-state, while the MAGIC detection corresponds to an exceptional, short-term flaring event during which no GeV -ray observatory was operating, but one may expect that the HE -ray flux at that time was larger than the Fermi-LAT high-state flux presented in Abdo et al. (2010) and shown in Figure 1. The radio emission from our model simulation is synchrotron-self-absorbed and therefore underpredicts the observed radio flux from 3C279. This suggests that the observed radio emission likely originates in more extended regions of the jet, beyond the radiative zone considered in our model.

In our simulation, the jet is — to within a factor of a few — in approximate equipartition between the powers carried in magnetic fields and in kinetic energy of particles: The power carried along the jet in the form of magnetic field (i.e., the Poynting Flux) is determined by

 LB=πR2Γ2βΓcB28π (35)

which, for our baseline fit to 3C279, yields  erg s. The particle kinetic luminosities in the observer’s frame are calculated from the equilibrium particle distributions as

 Li=πR2Γ2βΓcmic2∞∫1dγini(γi)γi (36)

where denotes the particle species considered. From numerically integrating the solution to the Fokker-Planck equation for both the proton and electron/positron distributions when equilibrium is reached, we find that the corresponding particle kinetic luminosities are  erg s and  erg s. With these values, the partition parameter between the combined particle kinetic luminosity and the power carried by the magnetic field, , where , is then . Our value for is similar to the values usually required by most previously published hadronic model interpretations of FSRQ SEDs. However, previously published works usually require parameters far out of equipartition. For example, in Böttcher et al. (2013), for their fit to the SED of 3C279, while our model produces a reasonable representation of the same SED with a jet near equipartition. This might be a consequence of the higher radiative efficiency in the parameter regime chosen here, with the inclusion of secondary muon and pion synchrotron radiation.

## 4 Simulated Lightcurves

We use the parameter set from our steady-state fit to the SED of 3C 279 as a baseline model from which we start out to apply perturbations to a few parameters in order to investigate the effects of these perturbations on the resulting multiwavelength light curves. Once the model has reached equilibrium as described in the previous section, we modify the time step to  s. This allows us to resolve light curve patterns on time scales characteristic for cooling effects of the relativistic protons. However, we are unable to diagnose the shorter-term variability, potentially caused by the radiative cooling of high energy electron-positron pairs generated from the decay of charged mesons, since their cooling time scales are significantly shorter than the size of the time step selected for these simulations. Note, again, that the implicit Crank-Nicholson scheme implemented for the solution of the Fokker-Planck equations guarantees that a stable solution for the electrons/positrons, muons, and pions is obtained even if the time step is longer than the radiative coolimg time scale. Decreasing the size of the time step would allow us to probe variability on those shorter time scales, but extending such simulations to time scales of the order of the proton cooling time scales would require prohibitively long simulation times.

After the simulation has reached equilibrium, one of the input parameters (, , , ) is modified in the form of a Gaussian perturbation in time. From the outputs produced in the simulations, we integrate the light curves in the following bandpasses: optical (R-band), X-ray (0.1 keV – 10 keV), HE -rays (20 MeV – 300 GeV) and VHE -rays (30 GeV – 100  TeV). The magnetic field perturbation is given by

 B(t)=B0+KB⋅e−(t−t0)2/2σ2 (37)

where  G denotes the equilibrium value for the magnetic field,  G parametrizes the amplitude of the perturbation, and and specify the time when the perturbation reaches its peak and the characteristic time scale of the perturbation, respectively. The chosen perturbation for the proton injection luminosity has the same functional form,

 Linj(t)=Linj,0+KL⋅e−(t−t0)2/2σ2 (38)

where  erg s is the equilibrium proton injection luminosity and is the amplitude of the perturbation. The perturbation of the acceleration time scale is chosen in such a way that the acceleration time scale decreases to a minimum during the peak of the perturbation. This is achieved with the following parametrization:

 tacc(t)=tacc,01+Kt⋅e−(t−t0)2/2σ2 (39)

where is the equilibrium value of the acceleration time scale and characterizes the amplitude of the perturbation. We also include a perturbation of the proton spectral index such that a flare is caused by a temporarily harder proton spectral index:

 qp(t)=qp,0−Kq⋅e−(t−t0)2/2σ2 (40)

where denotes the equilibrium value for the proton spectral index and denotes the strength of the perturbation. For all four perturbations, we choose a width of  s, corresponding to approximately 10 light-crossing time scales through the emission region, in the co-moving frame.

For the example of the modification of the magnetic field, Figure 2 illustrates snap-shot SEDs at various times throughout the simulation. The corresponding light curves (normalized to the respective peak fluxes) extracted from our time-dependent simulations are shown in figures 3 to 6.

A temporary increase in the magnetic field obviously leads to a marked increase in the proton synchrotron (primarily HE -rays) and electron-synchrotron (IR – optical) spectral components. The corresponding increase of the synchrotron-photon energy density leads to a larger pion-production (and subsequent pion- and muon-decay) rate. The resulting pions and muons are also subjected to the increased magnetic field, thus strongly increasing the contribution of muon and pion synchrotron to the SED. This increase in synchrotron emission from secondary particles leads to a distinct VHE -ray flare, slightly delayed with respect to the primary proton-synchrotron (HE -ray) flare. After the secondary particles have decayed to electrons and positrons, those secondaries also cool via synchrotron emission, producing a marked flare in the X-ray bandpass, with a very short (on the electron-synchrotron cooling time-scale) delay with respect to the VHE -ray flare. The enhanced pion and muon decay rates also lead to an increased neutrino flux, approximately coincident with the secondary electron/positron synchrotron (X-ray) flare.

A perturbation in the proton injection luminosity causes the proton synchrotron emission to increase producing primarily a HE -ray flare. This increase in both the number of protons and proton-synchrotron photons leads to a strongly enhanced pion (and subsequently, muon) production rate. The synchrotron emission from these additional high energy pions and muons produces a slightly delayed flare in the VHE -ray bandpass, in tandem with an increased neutrino flux from pion and muon decay. As in the case of the magnetic-field perturbation, the additional secondary electrons/positrons from the pion and muon decay then produce a delayed X-ray synchrotron flare. As these secondaries eventually cool down to even lower energies, their synchrotron emission contributes even to the optical (R-band) flux, leading to a slightly delayed flare in that band.

The perturbation characterized by an increase of the acceleration efficiency, leads to interesting features which are quite distinct from the B-field and injection-luminosity enhancements discussed above. With an increase in the stochastic acceleration efficiency, particles diffuse more efficiently to lower and higher energies. As protons are accelerated to higher energies, the proton synchrotron emission extends to higher energies, now entering the VHE -ray regime, leading to a prompt VHE -ray flare. The ultrarelativistic protons interact with the enhanced synchrotron radiation field, thus increasing the pion and muon production rates. The pions and muons themselves are subject to the increased acceleration efficiency and are thus pushed to higher energies, leading to a delayed, secondary VHE -ray flare due to pion and muon synchrotron radiation. All particle distributions cool due to synchrotron emission so that the spectral components gradually progress to lower frequencies. This leads to delayed flares in the HE -rays as well as X-rays and optical (R-band).

The perturbation of the proton spectral index also produces a interesting, distinct features. Due to the harder proton spectrum, the primary proton synchrotron emission temporarily also makes a larger contribution to the VHE -ray emission, leading to a brief, pronounced VHE flare. As in the case of the other perturbations discussed above, this also leads to increased pion and muon production rates, leading to delayed X-ray, optical, and neutrino flares. The synchrotron emission from the cooled, additional highest-energy protons produce a delayed flare in the HE -ray bandpass. As the flare progresses, the stronger proton cooling due to pion production results in a temporarily lower high-energy cut-off of the proton spectrum, because the high-energy cut-off of the proton injection spectrum remained unchanged, while radiative cooling becomes more efficient. As a result, the high-energy end of the proton synchrotron spectrum no longer makes a significant contribution to the VHE -ray flux, and the pion production rate (and subsequent pion and muon synchrotron emission) decreases temporarily. This leads to a dip in the VHE -ray light curve, before the perturbation subsides and radiative equilibrium is re-established.

The distinct features in these lightcurves can be used as a key diagnostic to differentiate between one-zone leptonic and hadronic models. In our previous study of the analogous flaring scenarios in a one-zone leptonic model in Diltz & Böttcher (2014), we predicted that a perturbation characterized by an increase in the electron acceleration efficiency would produce a deficit in the X-ray emission, while producing marked flares in the R-band, HE and VHE bandpasses. In leptonic models to the SEDs of FSRQs (such as 3C279), the X-rays are typically dominated by the low-energy end of the SSC emission. The drop in the X-ray flux was therefore attributed to a shift of the SSC emission to higher energies as a consequence of the increased electron acceleration efficiency. In contrast, in hadronic-model fits, the X-rays are dominated by synchrotron radiation of relativistic protons at intermediate energies (). When an acceleration-time-scale perturbation is applied to the hadronic model, all particle populations (including protons, pions, and muons) are accelerated to higher energies, without substantially affecting the particle distributions at intermediate energies. This causes increased pion, muon and electron-positron production rates. Following the subsequent radiative cooling of secondary electrons/positrons, their synchrotron emission leads to a delayed X-ray flare.

Additional marked differences are the delayed VHE -ray plateau found in our simulation of the acceleration-efficiency perturbation and the dip in the VHE light curve predicted for the proton spectral-index perturbation, both of which are not expected in leptonic models. These marked differences in the multiwavelength light curves may serve as diatnostics to distinguish between one-zone leptonic and hadronic models.

## 5 Discrete Correlation Analysis:

In order to further analyze cross-correlations and time lags between the simulated light curves discussed in the previous section, we calculate the discrete correlation functions (DCF, Edelson et al., 1988) between the light curves in the various bandpasses. In order to quantify the preferred values of the strength of the correlation (the maximum amplitude of the DCF) and inter-band time lag, a Gaussian fit to the DCFs was performed that minimized the chi-square between the data set and a fitting function of the form

 f(τ)=F1⋅e−(τ−τpk)2/2σ2 (41)

For this purpose, in order to be able to evaluate a value, we arbitrarily assumed a relative flux error of 1 % for each simulated light curve point when calculating the DCFs and their errors. In this discussion, we focus on the X-ray through -ray portion of the spectrum, and thus, on the DCFs between X-rays, HE -rays, and VHE -rays. This is largely motivated by significant differences between leptonic and hadronic models in the X-ray and -ray light curves for the acceleration time scale and proton spectral index perturbations. The best fit parameters for the various flaring scenarios are listed in Table 3.

The discrete correlation functions show strong correlations between the X-ray, HE, and VHE bandpasses for all perturbations considered in this paper. For both the magnetic field and proton injection luminosity perturbations, we find that the HE -ray flare is followed by a flare in VHE -rays and finally by a flare in the X-ray bandpass. This gives credence to the physical scenario discussed in the previous section that an increase in the synchrotron photon fields will subsequently generate flares in the VHE -ray bandpass and then a delayed X-ray flare.

For the acceleration timescale and the proton spectral index perturbation, the DCF analysis confirms the leading VHE -ray flare, followed by delayed HE -ray and X-ray flares. Time lags between the X-ray and -ray bands are typically – a few hours. Within error bars, the time lags determined from the DCFs agree with those extracted from visual inspection of the light curves. Variability on such time scales (and, thus, corresponding inter-band time lags) is easily measurable in X-rays and VHE -rays. However, the measurement of HE -ray variability on time scales of a few hours by Fermi-LAT may be possible only in extraordinarily high flux states.

## 6 Summary and Conclusions:

We have then simulated light curves by applying perturbations to a various input parameters in our code. The perturbations of the magnetic field and proton injection luminosity produced strong correlations between all bandpasses with hour time lags between the HE -ray and X-ray bandpass and hour time lags between the VHE and HE -ray bandpasses. Also a temporary increase in the stochastic acceleration efficiency leads to correlated flares in the -ray and X-ray bandpasses. This is in contrast to the the effects of such a perturbation on a time-dependent leptonic model (Diltz & Böttcher, 2014), in which a drop in X-ray emission was predicted. The predicted variability features are well within reach of observational capabilities of currently operating X-ray and VHE -ray observatories, but require extraordinarily high flux states to be measurable by Fermi-LAT. Our baseline (quiescent-state) model fit simulations predict integrated neutrino number fluxes at Earth, over the IceCube energy range for all three neutrino species, of  cm s. Given IceCube’s effective area of  cm, this predicts neutrino detection rates of  yr, thus requiring  years for a significant detection of neutrinos from 3C279 in quiescence. Even during flaring states, as studied in this paper, the neutrino flux is expected to increase by factors of a few – a few tens, to expected detection rates of  s, rendering it unlikely that IceCube would be able to detect a neutrino signal correlated with -ray flares from 3C279.

The most interesting features in our simulated lightcurve were plateaus and dips in the VHE -ray bandpass as a result of perturbations of the acceleration time scale or the proton injection spectral index. These plateaus are primarily caused by delayed synchrotron radiation from the secondary pions and muons. Such VHE light curve plateaus / dips are not predicted in one zone leptonic models and could be a tell tale signature of hadronic emission from blazar jets in parameter regimes in which muon and pion synchrotron emission is non-negligible.

## Acknowledgments

This work was funded by NASA through Astrophysics Data Analysis Program grant NNX12AE43G. The work of M.B. is supported through the South African Research Chair Initiative of the National Research Foundation and the Department of Science and Technology of South Africa, under SARChI Chair grant No. 64789.

## References

• Abdo et al. (2010) Abdo, A. A., et al., 2010, 716, 30
• Aharonian et al. (2000) Aharonian, F. A., 2000, New Astron., 5, 377
• Albert et al. (2008) Albert, J., et al., 2008, Science, 320, 1752
• Barr et al. (1988) Barr, S., Gaisser, T. K., et al., 1988, Phys. Letters 214, 147.
• Ballo et al. (2002) Ballo, L., Maraschi, L., et al., 2002, ApJ, 567, 50B
• Bentz et al. (2013) Bentz, M. C., Denney, K. D., et al., 2013, ApJ, 767, 149
• Blandford & Levinson (1995) Blandford, R. D. & Levinson, A., 1995, ApJ, 441, 79
• Blazejowski et al. (2000) Blazejowski, M., et al., 2000, ApJ, 545, 107
• Böttcher & Schlickeiser (1997) Böttcher, M., Schlickeiser, R., 1997, A&A, 325, 866
• Böttcher & Chiang (2002) Böttcher, M., Chiang, J., 2002, ApJ, 581, 127
• Böttcher et al. (2007) Böttcher, M., et al., 2007, ApJ, 670, 968
• Böttcher et al. (2009) Böttcher, M., Reimer, A., & Marscher, A. P., 2009, ApJ, 703, 1168
• Böttcher et al. (2012) Böttcher, M., Harris, D. E., Krawczynski, H., 2012, ”Relativistic Jets from Active Galactic Nuclei” WILEY-VCH Verlag GmBH & Co.
• Böttcher et al. (2013) Böttcher, M., Reimer, A., Sweeney, K., & Prakash, A., 2013, ApJ, 768, 54
• Dermer et al. (1992) Dermer, C. D., Schlickeiser, R., & Mastichiadis, A., 1992, A&A, 256, L27
• Dermer & Schlickeiser (1993) Dermer, C. D., Schlickeiser, R., 1993, ApJ, 416, 458
• Dermer & Menon (2009) Dermer, C. D., Menon, G., 2009, ”High Energy Radiation from Black Holes” Princeton University Press
• Diltz & Böttcher (2014) Diltz, C. & Böttcher, M., 2014, JHEAP, 63D
• Edelson et al. (1988) Edelson, R. A., Krolik, J. H., 1988, ApJ, 333, 646
• Hartman et al. (1996) Hartman, R. C., et al., 1996, ApJ, 461, 698
• Hartman et al. (2001) Hartman, R. C. et al., 2001, ApJ, 553, 683
• Hovatta et al. (2009) Hovatta, T., Valtaoja, E., Tornikoski, M., & Lähteenmk̈i, A., 2009, A&A, 498, 723
• Hümmer et al. (2010) Hümmer, S., Rüger, M., Spainer, F., Winter, W., 2010, ApJ 721 630H
• Jones (1968) Jones, F. C., 1968, Phys. Rev., 167, 1159
• Kelner & Aharonian (2008) Kelner, S. R., Aharonian, F. A., 2008, Phys. Rev. D78, 034013
• Mannheim & Biermann (1992) Mannheim. K., & Biermann, P. L., 1992, A&A, 253, L21
• Mannheim (1993) Mannheim, K., 1993, A&A, 269, 69
• Maraschi et al. (1994) Maraschi, L., et al., 1994, ApJ, 435, L91
• Mastichiadis & Kirk (1995) Mastichiadis, A., & Kirk. J. G., 1995, A&A, 295, 613
• Mastichiadis & Kirk (2005) Mastichiadis, A., & Kirk. J. G., 2005, A&A, 433, 765
• Mücke & Protheroe (2001) Mücke, A., & Protheroe, R. J., 2001, APh, 15, 121
• Mücke et al. (2003) Mücke, A., Protheroe, R. J., Engel, R., Rachen, J. P., & Stanev, T., 2003, APh, 18, 593
• Schlickeiser (1984a) Schlickeiser, R., 1984a, A&A, 136, 227
• Schlickeiser (2002) Schlickeiser, R., 2002, Astronomy & Astrophysics Library
• Sikora et al. (1994) Sikora, M., Begelman, M., & Rees, M., 1994 ApJ, 421, 153
• Urry (1998) Urry, C. M., 1998, Advances in Space Research, 21 ,89
• Weidinger & Spanier (2015) Weidinger, M., Spanier, F., 2015, A&A, 573, A7
• Woo et al. (2002) Woo, J. H., Urry, C. M., 2002, ApJ, 579, 530W
• Xie et al. (2008) Xie, Z. H., Hao, J. M., Du, L. M., Zhang, X., & Jia, Z. L., 2008, PSAP, 120, 477
• Zdziarski & Böttcher (2015) Zdziarski, A. A., & Böttcher, M., 2015, MNRAS, submitted (arXiv:1501.06124)
• Zhang & Böttcher (2013) Zhang, H., & Böttcher, M., 2013, ApJ, 774, 18
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters