TIME DEPENDENT HADRONIC MODELING OF FLAT SPECTRUM RADIO QUASARS

Blazars are a subcategory of active galactic nuclei (AGNs) 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 (HE) γ-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 infrared (IR) 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; DimitraKoudis et al. 2012; Petropoulou et al. 2014). 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 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 timescales expected for electrons/positrons versus 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 timescales 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, ${{\gamma }_{p,\,{\rm max} }}$, and magnetic field B to $B\;{{\gamma }_{p,\,{\rm max} }}\ll 5.6\times {{10}^{10}}$ G (Böttcher et al. 2013). If blazar jets are the sites of the acceleration of ultra-high-energy cosmic rays (${{E}_{p}}\gtrsim {{10}^{19}}$ eV), then such models are only applicable in the range of magnetic fields of $B\lesssim 5$ 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 light curves in the various bandpasses. We describe the model setup in Section 2; we present a model fit to the SED of the blazar 3C 279, using an equilibrium solution of our code, in Section 3; we then study the light curves resulting from the Gaussian perturbation of the selected parameters in Section 4, and compute the cross-correlation functions between the various light curves, in Section 5; we present a summary and brief discussion in Section 6.

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 R, 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, ${\Delta }{{t}_{\operatorname{var}}}$, through

$$\begin{eqnarray}&&R=\frac{c\cdot {\Delta }{{t}_{\operatorname{var}}}\cdot \delta }{1+z}\end{eqnarray} \tag{ 1 }$$

where z is the redshift to the source and δ is the Doppler factor. The time evolution of the proton distribution is described by a Fokker–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:

$$\begin{eqnarray}&&{{Q}_{p}}(\gamma )={{Q}_{p,0}}{{\gamma }^{-{{q}_{p}}}}H(\gamma ;{{\gamma }_{p,{\rm min} }},{{\gamma }_{p,{\rm max} }})\end{eqnarray} \tag{ 2 }$$

where $H(\gamma ;{{\gamma }_{p,{\rm min} }},{{\gamma }_{p,{\rm max} }})$ denotes the Heaviside function defined by H = 1 if ${{\gamma }_{p,{\rm min} }}\leqslant \gamma \leqslant {{\gamma }_{p,{\rm max} }}$ and H = 0 otherwise. The normalization factor for the proton injection spectrum is constrained through the proton injection luminosity, ${{L}_{p,{\rm inj}}}$:

$$\begin{eqnarray}{{Q}_{p,0}}=\left\{ \begin{array}{ccccccccccccccc} \frac{{{L}_{p,{\rm inj}}}}{{{V}_{b}}{{m}_{p}}{{c}^{2}}}\frac{2-q}{\gamma _{p,{\rm max} }^{2-q}-\gamma _{p,{\rm min} }^{2-q}} & {\rm if}\;q\ne 2, \\ \frac{{{L}_{{\rm p},\;{\rm inj}}}}{{{V}_{b}}{{m}_{p}}{{c}^{2}}{\rm ln} (\frac{{{\gamma }_{p,{\rm max} }}}{{{\gamma }_{p,{\rm min} }}})} & {\rm if}\;q=2. \\ \end{array} \right.\end{eqnarray} \tag{ 3 }$$

where V_b denotes the comoving volume of the emission region and m_p 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,

$$\begin{eqnarray}&&D(\gamma )=K\cdot {{\gamma }^{p}}\end{eqnarray} \tag{ 4 }$$

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

$$\begin{eqnarray}&&K=\frac{1}{(a+2){{t}_{{\rm acc}}}}\end{eqnarray} \tag{ 5 }$$

where $a=v_{s}^{2}/v_{A}^{2}$. We invoke a diffusion coefficient with a spectral index of p = 2 (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:

$$\begin{eqnarray}&&\begin{array}{l} {{Q}_{b}}({{E}_{b}})=\int _{{{E}_{b}}}^{\infty }\frac{d{{E}_{p}}}{{{E}_{p}}}{{N}_{p}}({{E}_{p}}) \\ \int _{\frac{{{\epsilon }_{{\rm th}}}{{m}_{p}}{{c}^{2}}}{2{{E}_{p}}}}^{\infty }d\epsilon {{n}_{\gamma }}(\epsilon ){{R}_{b}}(x,y) \\ \end{array}\end{eqnarray} \tag{ 6 }$$

where ${{N}_{p}}({{E}_{p}})$ denotes the proton distribution, ${{n}_{\gamma }}(\epsilon )$ denotes the photon field that the protons interact with as a function of the normalized photon energy $\epsilon =h\nu /({{m}_{e}}{{c}^{2}})$, and ${{\epsilon }_{{\rm th}}}=294$ (corresponding to an energy of 150 MeV) represents the threshold below which the cross sections are zero. The dimensionless variables x and y are given by

$$\begin{eqnarray}&&x=\frac{{{E}_{b}}}{{{E}_{p}}}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&y=\frac{{{E}_{p}}\epsilon }{{{m}_{p}}{{c}^{2}}}.\end{eqnarray} \tag{ 8 }$$

The response function $R(x,y)$ in the photo production rate of pions is given by

$$\begin{eqnarray}\begin{array}{rcl} {{R}_{b}}(x,y) & = & \mathop{\sum }\limits_{{\rm IT}}{{R}^{{\rm IT}}}(x,y) \\ {} & = & \mathop{\sum }\limits_{{\rm IT}}\frac{1}{2{{y}^{2}}}\int _{{{\epsilon }_{{\rm th}}}}^{2y}d{{\epsilon }_{r}}{{\epsilon }_{r}}{{\sigma }^{{\rm IT}}}({{\epsilon }_{r}}) \\ {} & {} & \times M_{b}^{{\rm IT}}({{\epsilon }_{r}})\delta (x-{{\chi }^{{\rm IT}}}({{\epsilon }_{r}})) \\ \end{array}\end{eqnarray} \tag{ 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, ${{\sigma }^{{\rm IT}}}({{\epsilon }_{r}})$. The functions $M_{b}^{{\rm IT}}({{\epsilon }_{r}})$ and ${{\chi }^{{\rm IT}}}({{\epsilon }_{r}})$ 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, ${{\epsilon }_{r}}$. The response function then simplifies to:

$$\begin{eqnarray}&&{{R}^{{\rm IT}}}(x,y)=\delta (x-{{\chi }^{{\rm IT}}})M_{b}^{{\rm IT}}{{f}^{{\rm IT}}}(y)\end{eqnarray} \tag{ 10 }$$

with

$$\begin{eqnarray}&&{{f}^{{\rm IT}}}(y)=\frac{1}{2{{y}^{2}}}\int _{{{\epsilon }_{{\rm th}}}}^{2y}d{{\epsilon }_{r}}{{\epsilon }_{r}}{{\sigma }^{{\rm IT}}}({{\epsilon }_{r}}).\end{eqnarray} \tag{ 11 }$$

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

$$\begin{eqnarray}&&\begin{array}{l} Q_{b}^{{\rm IT}}={{N}_{p}}(\frac{{{E}_{b}}}{{{\chi }^{{\rm IT}}}})\cdot \frac{{{m}_{p}}{{c}^{2}}}{{{E}_{b}}} \\ \quad \int _{{{\epsilon }_{{\rm th}}}/2}^{\infty }dy{{n}_{\gamma }}(\frac{{{m}_{p}}{{c}^{2}}\;y{{\chi }^{{\rm IT}}}}{{{E}_{b}}})M_{b}^{{\rm IT}}{{f}^{{\rm IT}}}(y). \\ \end{array}\end{eqnarray} \tag{ 12 }$$

This single integral is easy to evaluate numerically. The photo-production rate of pions now depends on the response function, ${{f}^{{\rm IT}}}(y)$, the multiplicities, M_b^IT, and the mean energy fraction deposited into the secondary particles, ${{\chi }^{{\rm IT}}}$, 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:

$$\begin{eqnarray}&&t_{{\rm cool}}^{-1}({{E}_{p}})=\mathop{\sum }\limits_{{\rm IT}}M_{p}^{{\rm IT}}{{{\Gamma }}^{{\rm IT}}}({{E}_{p}}){{K}^{{\rm IT}}}\end{eqnarray} \tag{ 13 }$$

where K^IT is the inelasticity and ${{{\Gamma }}^{{\rm IT}}}({{E}_{p}})$ the interaction rate of a given interaction type, given by:

$$\begin{eqnarray}&&{{{\Gamma }}^{{\rm IT}}}({{E}_{p}})=\int _{\frac{{{\epsilon }_{{\rm th}}}{{m}_{p}}{{c}^{2}}}{2{{E}_{p}}}}^{\infty }d\epsilon {{n}_{\gamma }}(\epsilon ){{f}^{{\rm IT}}}(\frac{{{E}_{p}}\epsilon }{{{m}_{p}}{{c}^{2}}}).\end{eqnarray} \tag{ 14 }$$

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

$$\begin{eqnarray}&&t_{{\rm esc},n}^{-1}({{E}_{p}})=\mathop{\sum }\limits_{{\rm IT},{{p}^{\prime }}\ne p}M_{{{p}^{\prime }}}^{{\rm IT}}{{{\Gamma }}^{{\rm IT}}}({{E}_{p}})\end{eqnarray} \tag{ 15 }$$

where p' 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:

$$\begin{eqnarray}&&{{\dot{\gamma }}_{p\gamma }}=-{{E}_{p}}\cdot t_{{\rm cool}}^{-1}({{E}_{p}})=-{{E}_{p}}\cdot \mathop{\sum }\limits_{{\rm IT}}M_{p}^{{\rm IT}}{{{\Gamma }}^{{\rm IT}}}({{E}_{p}}){{K}^{{\rm IT}}}.\end{eqnarray} \tag{ 16 }$$

Given that the radiative cooling timescales 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 $\theta \sim 1/{\Gamma }$, the adiabatic cooling rate is ${{\dot{\gamma }}_{{\rm ad}}}=-3c\theta \gamma /R$. 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):

$$\begin{eqnarray}&&\begin{array}{l} \frac{\partial {{n}_{p}}(\gamma ,t)}{\partial t}=\frac{\partial }{\partial \gamma }[\frac{{{\gamma }^{2}}}{(a+2){{t}_{{\rm acc}}}}\frac{\partial {{n}_{p}}(\gamma ,t)}{\partial \gamma }] \\ \quad -\frac{\partial }{\partial \gamma }(\dot{\gamma }\cdot {{n}_{p}}(\gamma ,t))+{{Q}_{p}}(\gamma ,t) \\ \quad -\frac{{{n}_{p}}(\gamma ,t)}{{{t}_{{\rm esc}}}}-\frac{{{n}_{p}}(\gamma ,t)}{{{t}_{{\rm esc},\;{\rm n}}}} \\ \end{array}\end{eqnarray} \tag{ 17 }$$

where ${{t}_{{\rm esc}}}$ denotes the dynamical escape time scale for the protons which we parameterize as a multiple of the light crossing time: ${{t}_{{\rm esc}}}=\eta R/c$ where $\eta \geqslant 1$. The value of η is kept as a free parameter. The term $\dot{\gamma }$ denotes the combined loss rates on the proton distribution due to adiabatic, synchrotron and pion production processes. The synchrotron loss rate is given by

$$\begin{eqnarray}&&{{\dot{\gamma }}_{{\rm rad}}}=-\frac{4}{3}c{{\sigma }_{T}}\frac{{{u}_{B}}}{{{m}_{e}}{{c}^{2}}}{{(\frac{{{m}_{e}}}{{{m}_{p}}})}^{3}}{{\gamma }^{2}}\end{eqnarray} \tag{ 18 }$$

where u_B 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 $t_{{\rm decay}}^{\prime }=2.8\times {{10}^{-17}}$ 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 ($t_{{\rm decay}}^{\prime }=2.6\times {{10}^{-8}}$ 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:

$$\begin{eqnarray}&&\begin{array}{l} \frac{\partial {{n}_{\pi }}(\gamma ,t)}{\partial t}=\frac{\partial }{\partial \gamma }[\frac{{{\gamma }^{2}}}{(a+2){{t}_{{\rm acc}}}}\frac{\partial {{n}_{\pi }}(\gamma ,t)}{\partial \gamma }] \\ \quad -\frac{\partial }{\partial \gamma }({{{\dot{\gamma }}}_{{\rm rad}}}{{n}_{p}}(\gamma ,t))+{{Q}_{\pi }}(\gamma ,t) \\ \quad -\frac{{{n}_{\pi }}(\gamma ,t)}{{{t}_{{\rm esc}}}}-\frac{{{n}_{\pi }}(\gamma ,t)}{\gamma t_{{\rm decay}}^{\prime }}. \\ \end{array}\end{eqnarray} \tag{ 19 }$$

The only major difference comes from the loss term due to the decay of charged pions with a characteristic timescale ${{t}_{{\rm decay}}}=\gamma t_{{\rm decay}}^{\prime }$. 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:

$$\begin{eqnarray}&&{{\pi }^{+}}\to {{\mu }^{+}}+{{\nu }_{\mu }}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{{\pi }^{-}}\to {{\mu }^{-}}+{{\bar{\nu }}_{\mu }}.\end{eqnarray} \tag{ 21 }$$

The decay term in the pion Fokker–Planck equation (last term in Equation (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:

$$\begin{eqnarray}&&{{\mu }^{+}}\to {{e}^{+}}+{{\nu }_{e}}+{{\bar{\nu }}_{\mu }}\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{{\mu }^{-}}\to {{e}^{-}}+{{\bar{\nu }}_{e}}+{{\nu }_{\mu }}.\end{eqnarray} \tag{ 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 $\gamma \gamma$ absorption and pair production, which allows us to properly follow the evolution of ultra-high-energy γ-ray induced pair cascades (see, e.g., Böttcher et al. 2013).

2.4. Radiative Contributions

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 ${{n}_{i}}(\gamma )$ of charged particles i within a tangled magnetic field is given by:

$$\begin{eqnarray}&&{{j}_{i,{\rm syn}}}(\nu ,t)=\frac{1}{4\pi }\int _{0}^{\infty }d\gamma {{n}_{i}}(\gamma ,t)\cdot {{P}_{i,\nu }}(\gamma )\end{eqnarray} \tag{ 24 }$$

where the term ${{P}_{i,\nu }}(\gamma )$ denotes the synchrotron power per unit frequency produced by a single charged particle of species i, and can be well approximated by (Böttcher et al. 2012):

$$\begin{eqnarray}&&\begin{array}{l} {{P}_{i,\nu }}(\gamma )=\frac{32\pi c}{9{\Gamma }(4/3)}\;r_{e}^{2}{{(\frac{{{m}_{e}}}{{{m}_{i}}})}^{2}} \\ \quad \cdot {{u}_{B}}\;{{\gamma }^{2}}\;\frac{{{\nu }^{1/3}}}{\nu _{c}^{4/3}}{{e}^{-\nu /{{\nu }_{c}}}} \\ \end{array}\end{eqnarray} \tag{ 25 }$$

where u_B denotes the energy density of the magnetic field, r_e the classical electron radius, and m_i the mass of a particle of species i. The critical synchrotron frequency ${{\nu }_{c}}$, is given by

$$\begin{eqnarray}&&{{\nu }_{c}}=4.2\times {{10}^{6}}\cdot B({\rm G})\cdot ({{m}_{e}}/{{m}_{i}})\cdot {{\gamma }^{2}}\;{\rm Hz}.\end{eqnarray} \tag{ 26 }$$

The synchrotron spectrum represents one component of the combined photon field that also includes the SSC and EC 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):

$$\begin{eqnarray}&&\begin{array}{l} \frac{\partial {{n}_{{\rm ph}}}(\nu ,t)}{\partial t}=\frac{4\pi }{h\nu }\cdot \mathop{\sum }\limits_{k}{{j}_{k,\nu }}(t) \\ \quad -{{n}_{{\rm ph}}}(\nu ,t)\cdot (\frac{1}{{{t}_{{\rm esc},\;{\rm ph}}}}+\frac{1}{{{t}_{{\rm abs}}}}) \\ \end{array}\end{eqnarray} \tag{ 27 }$$

where the sum is over the all radiation mechanisms and all particle species, ${{t}_{{\rm esc},\;{\rm ph}}}=4R/3c$ denotes the photon escape time scale, and ${{t}_{{\rm abs}}}$ denotes the absorption time scale due to synchrotron-self-absorption by electrons and $\gamma \gamma$ absorption. The absorption time scale can be defined through the opacity as

$$\begin{eqnarray}&&{{t}_{{\rm abs}}}=\frac{R}{c\cdot ({{\tau }_{{\rm SSA}}}+{{\tau }_{\gamma \gamma }})}\end{eqnarray} \tag{ 28 }$$

where ${{\tau }_{{\rm SSA}}}$ and ${{\tau }_{\gamma \gamma }}$ denote the synchrotron-self-absorption and $\gamma \gamma$ 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 $\gamma \gamma$ 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:

$$\begin{eqnarray}&&\nu F_{\nu }^{{\rm obs}}({{\nu }^{{\rm obs}}},{{t}^{{\rm obs}}})=\frac{h\cdot {{\nu }^{2}}\cdot {{n}_{{\rm ph}}}(\nu ,t)\cdot {{\delta }^{4}}\cdot {{V}_{b}}}{4\pi d_{L}^{2}\cdot {{t}_{{\rm esc},{\rm ph}}}}\end{eqnarray} \tag{ 29 }$$

with ${{\nu }^{{\rm obs}}}=\delta \nu$ and ${\Delta }{{t}^{{\rm obs}}}={\Delta }t/\delta$.

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, ${{D}_{\pi }}({{E}_{\pi }})$, which is given by the decay term in the pion Fokker–Planck equation:

$$\begin{eqnarray}&&{{D}_{\pi }}({{E}_{\pi }})=\frac{{{n}_{\pi }}(\gamma ,t)}{t_{{\rm decay}}^{\prime }\;\gamma }.\end{eqnarray} \tag{ 30 }$$

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

$$\begin{eqnarray}&&{{Q}_{\nu }}({{E}_{\nu }})=\int _{{{E}_{\nu }}{{(1-{{r}_{M}})}^{-1}}}^{\infty }\frac{d{{E}_{\pi }}}{{{E}_{\pi }}}\cdot \frac{{{D}_{\pi }}({{E}_{\pi }})}{1-{{r}_{M}}}\end{eqnarray} \tag{ 31 }$$

where ${{r}_{M}}=m_{\mu }^{2}/m_{\pi }^{2}$.

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:

$$\begin{eqnarray}&&{{Q}_{\nu }}({{E}_{\nu }})=\int _{{{E}_{\nu }}}^{\infty }d{{E}_{\mu }}{{D}_{\mu }}({{E}_{\mu }})\frac{dn}{d{{E}_{\nu }}}.\end{eqnarray} \tag{ 32 }$$

Using the dimensionless scalar variable $m={{E}_{\nu }}/{{E}_{\mu }}$, we can recast Equation (32) into the form:

$$\begin{eqnarray}&&{{Q}_{\nu }}({{E}_{\nu }})=\int _{0}^{1}dm\frac{{{D}_{\mu }}({{E}_{\nu }}/m)}{m}\cdot \frac{dn}{dm}\end{eqnarray} \tag{ 33 }$$

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

$$\begin{eqnarray}&&\frac{dn}{dm}={{g}_{0}}(m)+{{g}_{1}}(m).\end{eqnarray} \tag{ 34 }$$

The scalar functions ${{g}_{0}}(m)$ and ${{g}_{1}}(m)$ 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.

Table 1.  Laboratory-frame Electron and Muon Neutrino Distribution Functions (from Barr et al. 1988)

	${{g}_{0}}(m)$	${{g}_{1}}(m)$
${{\nu }_{\mu }}\;:$	$\frac{5}{3}-3{{m}^{2}}+\frac{4}{3}{{m}^{3}}$	$\frac{1}{3}-3{{m}^{2}}+\frac{8}{3}{{m}^{3}}$
${{\nu }_{e}}\;:$	$2-6{{m}^{2}}+4{{m}^{3}}$	$-2+12m-18{{m}^{2}}+8{{m}^{3}}$

Download table as:  ASCIITypeset image

3C 279 $(z=0.538)$ 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 3C 279 (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 10⁷ s. This time step size is considerably longer than the timescales 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 3C 279, we have the following observational parameters (see Böttcher et al. 2013 for references to the observational data): z = 0.536, ${{\beta }_{\bot ,{\rm app}}}=20.1$ (the apparent transverse velocity of individual jet components, normalized to the speed of light), ${\Delta }{{t}_{\operatorname{var}}}\sim 2$ days, ${{L}_{{\rm disk}}}=2.0\times {{10}^{45}}$ erg s⁻¹, and ${{L}_{{\rm BLR}}}=2.0\times {{10}^{44}}$ erg s⁻¹. The superluminal motion speed sets a lower limit to the bulk Lorentz factor, ${\Gamma }\gt 21$. The observing angle is set by using the relation ${{\theta }_{{\rm obs}}}=1/{\Gamma }$ so that $\delta ={\Gamma }$. From the variability time scale, we can constrain the location of the emission region along the jet, ${{R}_{{\rm axis}}}\sim 2\;{{{\Gamma }}^{2}}\;c\;{{t}_{v}}/(1+z)\approx {{10}^{18}}\;{\rm cm}$. 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 (4–8) $\times \;{{10}^{8}}\;{{M}_{{\rm sol}}}$ (Woo & Urry 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 $5.0\times {{10}^{3}}\;K$, 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.

Table 2.  Parameter Values Used for the Equilibrium Fit to the SED of 3C 279

Parameter	Symbol	Value
Magnetic field	B	150 G
Radius of emission region	R	8.69 × 1015 cm
Constant multiple for escape time scale	η	6.0
Bulk Lorentz factor	${\Gamma }$	21
Observing angle	θobs	4.76 × 10−2 rad
Minimum proton Lorentz factor	γp,min	1.0
Maximum proton Lorentz factor	γp,max	4.5 × 108
Proton injection spectral index	qp	2.2
Proton injection luminosity	Lp,inj	$3.5\;\times \;{{10}^{46}}$ erg s−1
Minimum electron Lorentz factor	γe,min	5.1 × 102
Maximum electron Lorentz factor	γe,max	1.0 × 104
Electron injection spectral index	qe	3.2
Electron injection luminosity	Le,inj	7.8 × 1041 erg s−1
Supermassive black hole mass	MBH	6.0 × 108 M⊙
Eddington ratio	lEdd	1.18 × 10−2
Blob location along the jet axis	Raxis	0.279 pc
Radius of the BLR	R'ext	0.071 pc
Energy density of the BLR in comoving frame	u'ext	3.68 × 10−4 erg cm−3
Blackbody temperature of BLR	TBB	5000 K
Ratio between the acceleration and escape time scales	tacc/tesc	32.5

Download table as:  ASCIITypeset image

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 $B{{\gamma }_{{\rm p},\;{\rm max} }}\gtrsim 5\times {{10}^{10}}\;{\rm G}$, 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 3C 279. This suggests that the observed radio emission likely originates in more extended regions of the jet, beyond the radiative zone considered in our model.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{{L}_{B}}=\pi {{R}^{2}}{{{\Gamma }}^{2}}{{\beta }_{{\Gamma }}}c\frac{{{B}^{2}}}{8\pi }\end{eqnarray} \tag{ 35 }$$

which, for our baseline fit to 3C 279, yields ${{L}_{B}}=2.8\times {{10}^{48}}$ erg s⁻¹. The particle kinetic luminosities in the observer's frame are calculated from the equilibrium particle distributions ${{n}_{i}}({{\gamma }_{i}})$ as

$$\begin{eqnarray}&&{{L}_{i}}=\pi {{R}^{2}}{{{\Gamma }}^{2}}{{\beta }_{{\Gamma }}}c\;{{m}_{i}}{{c}^{2}}\int _{1}^{\infty }d{{\gamma }_{i}}\;{{n}_{i}}({{\gamma }_{i}})\;{{\gamma }_{i}}\end{eqnarray} \tag{ 36 }$$

where i 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 ${{L}_{p}}=9.7\times {{10}^{48}}$ erg s⁻¹ and ${{L}_{e}}=3.5\times {{10}^{43}}$ erg s⁻¹. With these values, the partition parameter between the combined particle kinetic luminosity and the power carried by the magnetic field, ${{\epsilon }_{B}}\equiv {{L}_{B}}/{{L}_{{\rm kin}}}$, where ${{L}_{{\rm kin}}}={{L}_{e}}+{{L}_{p}}$, is then ${{\epsilon }_{B}}\approx 0.29$. Our value for L_p 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), ${{L}_{B}}/{{L}_{p}}=7.9\times {{10}^{-3}}$ for their fit to the SED of 3C 279, 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.

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 ${\Delta }t=2.0\times {{10}^{4}}$ s. This allows us to resolve light curve patterns on timescales 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 timescales 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 timescales, but extending such simulations to timescales of the order of the proton cooling timescales would require prohibitively long simulation times.

After the simulation has reached equilibrium, one of the input parameters (B, ${{L}_{{\rm p},\;{\rm inj}}}$, ${{t}_{{\rm acc}}}$, ${{q}_{{\rm p}}}$) 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–10 keV), HE γ-rays (20 MeV–300 GeV), and VHE γ-rays (30 GeV–100 TeV). The magnetic field perturbation is given by

$$\begin{eqnarray}&&B(t)={{B}_{0}}+{{K}_{B}}\cdot {{e}^{-{{(t-{{t}_{0}})}^{2}}/2{{\sigma }^{2}}}}\end{eqnarray} \tag{ 37 }$$

where ${{B}_{0}}=150$ G denotes the equilibrium value for the magnetic field, ${{K}_{B}}=250$ G parametrizes the amplitude of the perturbation, and t₀ 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,

$$\begin{eqnarray}&&{{L}_{{\rm inj}}}(t)={{L}_{{\rm inj},\;0}}+{{K}_{L}}\cdot {{e}^{-{{(t-{{t}_{0}})}^{2}}/2{{\sigma }^{2}}}}\end{eqnarray} \tag{ 38 }$$

where ${{L}_{{\rm inj},\;0}}=3.5\times {{10}^{46}}$ erg s⁻¹ is the equilibrium proton injection luminosity and ${{K}_{L}}=0.3\times {{L}_{{\rm inj},\;0}}$ 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:

$$\begin{eqnarray}&&{{t}_{{\rm acc}}}(t)=\frac{{{t}_{{\rm acc},\;0}}}{1+{{K}_{t}}\cdot {{e}^{-{{(t-{{t}_{0}})}^{2}}/2{{\sigma }^{2}}}}}\end{eqnarray} \tag{ 39 }$$

where ${{t}_{{\rm acc},\;0}}$ is the equilibrium value of the acceleration time scale and ${{K}_{t}}=14$ 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:

$$\begin{eqnarray}&&{{q}_{p}}(t)={{q}_{p,0}}-{{K}_{q}}\cdot {{e}^{-{{(t-{{t}_{0}})}^{2}}/2{{\sigma }^{2}}}}\end{eqnarray} \tag{ 40 }$$

where ${{q}_{p,0}}$ denotes the equilibrium value for the proton spectral index and ${{K}_{q}}=1.0$ denotes the strength of the perturbation. For all four perturbations, we choose a width of $\sigma =1.0\times {{10}^{5}}$ s, corresponding to approximately 10 light-crossing timescales 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–6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The distinct features in these light curves 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 3C 279), 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 (${{\gamma }_{p}}\sim {{10}^{6}}$). 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.

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 & Krolik 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

$$\begin{eqnarray}&&f(\tau )={{F}_{1}}\cdot {{e}^{-{{(\tau -{{\tau }_{{\rm pk}}})}^{2}}/2{{\sigma }^{2}}}}.\end{eqnarray} \tag{ 41 }$$

For this purpose, in order to be able to evaluate a ${{\chi }^{2}}$ 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.

Table 3.  Best-fit DCF Correlation Strengths and Time Lags from the Gaussian Fits to the Discrete Correlation Functions

Bands	Scenario	F1	$\sigma ({\rm s})$	${{\tau }_{{\rm pk}}}({\rm s})$	Figures
X-HE	B	1.01	(1.59±0.13)×104	(−9.88±1.03)×103	7
HE-VHE	B	0.98	(1.67±0.19)×104	(−7.32±1.29)×103	8
X-HE	Lp,inj	1.04	(1.69±0.15)×104	(−1.14±0.11)×104	9
HE-VHE	Lp,inj	1.01	(1.69±0.21)×104	(−3.19±1.34)×103	10
X-HE	tacc	0.99	(1.49±0.17)×104	(1.83±1.19)×103	11
HE-VHE	tacc	0.91	(1.42±0.18)×104	(1.46±0.14)×104	12
X-HE	qp	0.99	(3.39±0.63)×104	(−4.08±2.31)×103	13
HE-VHE	qp	0.87	(3.43±1.05)×104	(2.64±0.53)×103	14

Download table as:  ASCIITypeset image

The DCFs 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 ∼1—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 timescales (and, thus, corresponding inter-band time lags) is easily measurable in X-rays and VHE γ-rays. However, the measurement of HE γ-ray variability on timescales of a few hours by Fermi-LAT may be possible only in extraordinarily high flux states.

In this paper, we have described a new time-dependent lepto-hadronic model for the broadband emission of relativistic jet sources (especially blazars) that incorporates the synchrotron radiation from secondary pions and muons, generated in photo-hadronic interactions. We use an equilibrium solution of our model to produce a rough fit-by-eye to the high-state SED of the FSRQ 3C 279. The broadband emission from infrared to VHE γ-rays can be explained by a combination of synchrotron emission from electrons, protons, pions, and muons. The parameters used for the fit are chosen so that the radiative contributions from the pion and muon synchrotron radiation are non-negligible, and their contribution is essential for an adequate fit of the unusually hard VHE γ-ray spectrum measured by MAGIC. Our fits for 3C 279 can be achieved with the jet being close to equipartition between the power carried in magnetic fields (Poynting flux) and the kinetic energy in protons and electrons (${{\epsilon }_{B}}\equiv {{L}_{B}}/{{L}_{{\rm kin}}}=0.29$). This contrasts most other hadronic models in which the particle kinetic luminosity is a few orders of magnitude larger than the magnetic luminosity. However, our model requires similarly large jet powers as other published hadronic blazar model fits (e.g., Böttcher et al. 2013, and references therein), greatly in excess of the Eddington luminosity of the central black hole in 3C 279 (${{L}_{{\rm Edd}}}\sim 8\times {{10}^{46}}$ erg s⁻¹). Nevertheless, as the jet power is strongly beamed perpendicular to the accretion flow, it does not provide a radiation pressure that would be able to shut off the accretion flow. Therefore, the Eddington limit argument may not apply in such a case. Nevertheless, the extreme jet production efficiency required for hadronic blazar jet models in general, may constitute a problem for this kind of models: the total jet power of ${{L}_{j}}\sim 1.3\times {{10}^{49}}$ erg s⁻¹ exceeds the radiative luminosity of the accretion disk (${{L}_{d}}\sim 2\times {{10}^{45}}$ erg s⁻¹) by almost 4 orders of magnitude. No jet production mechanism is currently known that would be able to produce steady jets with this efficiency; however, the current understanding of the production of relativistic jets from supermassive black holes is still very limited. This issue is further discussed in Zdziarski & Böttcher (2015).

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 3–4 hr time lags between the HE γ-ray and X-ray bandpass and 1–2 hr 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 $\approx {{10}^{-16}}$ cm⁻² s⁻¹. Given IceCube's effective area of ${{A}_{{\rm eff}}}(\gt 100\;{\rm TeV})\approx {{10}^{8}}$ cm²,this predicts neutrino detection rates of ∼0.3 yr⁻¹, thus requiring $\gtrsim 10$ yr for a significant detection of neutrinos from 3C 279 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 $\sim {{10}^{-7}}$ s⁻¹, rendering it unlikely that IceCube would be able to detect a neutrino signal correlated with γ-ray flares from 3C 279.

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.

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.