© ESO, 2014

1. Introduction

The origin of the characteristic highly non-thermal radiation from blazars, a subclass of active galactic nuclei (AGN) with the highly relativistic outflow emerging under a small angle to the line of sight, is undoubtedly the jet. The two prominent humps in the spectral energy distribution (SED) of these objects, one in the optical to X-ray regime and a second one occurring at the highest energies up to TeV, require high beaming factors which only the jet can provide (Urry 1998). Blazars have gained a lot of interest over the past decades mainly as a result of the discovery of their variable very high energy (VHE) emission with variations of many orders of magnitudes in the TeV regime using Air-Cherenkov telescopes like H.E.S.S., MAGIC, and VERITAS entrenching gamma-ray astronomy as its own field of research, from Punch et al. (1992) to Aharonian et al. (2007). With the Fermi satellite and its all-sky survey capabilities (in orbit since 2009), the availability of multiwavelength (MWL) data (combined with radio, optical, and X-ray observations using the Swift satellite, RXTE, Chandra, etc.) has rapidly increased so that the variable emission of AGN can be monitored in different energy bands. This data and the steadily increasing number of discovered blazars provides a solid basis for systematic investigations with theoretical emission models.

There have been many attempts to classify blazars, for example using visible line emission (Strittmatter et al. 1972; Marcha et al. 1996; Urry 1998), which in the end was summarized in the blazar sequence, first proposed by Fossati et al. (1998). Although there are exceptions, the most luminous blazars, flat spectrum radio quasars (FSRQ), peak at the lowest energies whereas the faintest, high frequency peaked BL Lac objects (HBL), show the highest peak frequencies. In between one finds low- and intermediate peaked BL Lac objects (LBL and IBL, respectively).

Recent modelling has shown that the SEDs from HBLs to IBLs are reliably reproduced by simple synchrotron self-Compton emission, i.e. the first peak in the spectrum is essentially synchrotron radiation from relativistic primary electrons and the second one arises from Compton upscattering of this synchrotron radiation by the very same electrons (Böttcher et al. 1997; Kirk et al. 1998; Tavecchio et al. 2001; Weidinger et al. 2010), but especially, when FSRQs are considered the origin of the second hump in the SED is still under debate, whether it is due to comptonized radiation from external broad line region photons, torus photons, or even accretion disk photons entering the jet (Dermer & Schlickeiser 1993; Böttcher & Dermer 1998) or due to highly relativistic primary protons within the jet and their synchrotron radiation as well as emerging photo-hadronic cascades (Mannheim 1993; Mücke et al. 2003; Böttcher et al. 2009). In addition to the stability and formation or the structure of the jet, which cannot be addressed with emission models, this composition of the jet and its difference in FSRQs and HBLs is one of the major open questions concerning blazars, and may provide a natural physical explanation of the phenomenological blazar sequence and its recently appearing envelope structure (Meyer et al. 2011).

In this paper we present a time-dependent and fully self-consistent, hybrid emission model for blazars where low energetic electrons and protons are co-accelerated via Fermi-I and Fermi-II processes within a confined region to high energies. Electrons lose energy mainly to the synchrotron and inverse Compton channels, protons radiate synchrotron photons and are subject to photo-meson production, with γγ-pair-production of highly energetic photons coupling the equations in a non-linear way. Proton collisions of the non-thermal proton distribution are irrelevant at typical densities within the jet and are neglected since no thermal background protons as an additional parameter (e.g. Eichmann et al. 2012) shall be considered. There have been very recent observations by Agudo et al. (2011) constraining the γ-ray emitting site of the jet a few pc away from the black hole. This strongly favours an emission scenario independent of external sources. We are able to model HBLs in the SSC case by setting the number density of injected protons Q_0,p⁺ → 0 (Weidinger et al. 2010; Weidinger & Spanier 2010) as well as FSRQs with Q_0,p⁺ ≠ 0 (Spanier & Weidinger 2011) with the very same model. Since our model is time-dependent we are able to exploit outbursts of blazars and the timing signatures in different energy bands to narrow down the parameters used in the modelling process and distinguish between leptonic and hadronic dominated jets.

In the next section we give a description of our model and its principle properties and assumptions, before we apply it to the recently detected HBL/IBL 1 ES 1011+496 which, already in the steady state, favours a very high magnetic field being present within the jet making it an outstanding source. In section 4 other, more general, physical implications are discussed.

Fig. 1 Basic model geometry. All escaping particles of the acceleration zone serve as injection for the radiation zone, i = e± for electrons/positrons and i = p+ for protons.

2. The model

The general modelsetup is assumed to be spherical with an acceleration zone nested inside a radiation zone. As the emitting region (blob) moves down the jet axis towards the observer with a bulk Lorentzfactor Γ, upstream material (electrons and protons) is picked up and the highly turbulent acceleration zone forms at the edge of the blob. Here both particle species undergo Fermi-I and Fermi-II processes up to relativistic energies with synchrotron losses in a turbulent magnetic field balancing the acceleration with regard to the maximum energy. On the far side of the blob, namely in the considerably larger radiation zone, acceleration is assumed to be inefficient. All calculations are conveniently made in the rest-frame of the blob, the geometry is shown in Fig. 1. Both zones are assumed to be homogeneous and to contain isotropic particle distributions as well as a randomly orientated magnetic field B, not to be confused with the helical magnetic field stabilizing the jet against the kinetic pressure of the plasma. The kinetic equations in each zone, one for each particle species i, are derived from the relativistic Vlasov equation (Schlickeiser 2002) applying the one-dimensional diffusion approximation in the highly relativistic p_i = γ_im_ic case. Furthermore, the hard-sphere approximation for the spatial diffusion coefficient occurring in the equations is used, see Weidinger et al. (2010). We note that when only electrons are picked up into the acceleration zone (Q_p⁺ → 0) the model reduces to the SSC case described in Weidinger et al. (2010) and Weidinger & Spanier (2010), which is an extension to Kirk et al. (1998). The acceleration timescales can be translated into the microphysics of the jet, being proportional to the particle’s mass in the energy independent case (Spanier & Weidinger 2011) (1)with the parallel spatial diffusion coefficients K_{| | ,i}, and v_s and v_A as shock and Alfvén speeds, respectively, providing the scattering centres mandatory for diffusive shock acceleration (DSA). Hence, Eq. (1)can be used to cross check the parameter typical values of the jet’s microphysics. Motivated by the typical gyro-timescale, the acceleration timescale for one species is set to be constant and proportional to the particle’s mass, while the timescale for the second species results naturally from Eq. (1). Following Vukcevic & Schlickeiser (2007) this ensures that the particles are accellerated well within the acceleration site. We note that this simple energy independent assumption may slightly underestimate the acceleration efficiency at small values of γ_i, but guarantees that energy gains at the highest values are not overestimated in terms of DSA. Eventually all particles escaping the acceleration zone enter the radiation zone downstream the jet. To ensure power-law spectra, as expected from shock acceleration, the escape timescale for the acceleration site is set to be constant and within the order of the acceleration timescale t_esc,i ∝ t_acc,i (Kirk et al. 1998; Weidinger et al. 2010). Because of the strong confinement of the particles in the radiation zone, as described in other models, one needs to account for all possible radiation mechanisms, not only the dominating synchrotron losses, as in the acceleration zone. This requires magnetic fields of G), typical for hadronic models (Mannheim 1993; Böttcher et al. 2009), whereas leptonic models typically assume G) (Weidinger et al. 2010; Tavecchio et al. 2001). Thus, the same order of equipartition, typically of %) in our model, regardless which regime (hybrid or leptonic only) considered, is reached. Hence, self-consistency commands that magnetic fields in blazar models cannot be arbitrarily high without considering non-thermal protons and their radiative output.

Because of the model geometry, the acceleration zone is assumed not to contribute to the model SED directly, hence we only solve the kinetic equations for the primary particles. The two relevant kinetic equations for the particles’ energy and volume density, which in the isotropic diffusive case take the form (2)with i being i = e⁻ and i = p⁺ for primary electrons and protons, respectively, and the synchrotron . Consequently a ∝ v_s/v_A is the ratio of Shock speed to Alfvén velocity, hence denoting the dominance of Fermi-I over Fermi-II processes in our model. Monoenergetic particles are injected into the radiation zone (3)to model the cumulation of particles with densities Q_0,i from the upstream direction with Lorentzfactors of γ_0,i by the blob. The density evolution in the acceleration zone starting with n_i = 0 at t₀ = 0 can be found in Fig. 2; the parameters are B = 10 G, , t_acc,e = 3 × 10³ s and a → ∞. As one would expect the slope in Fig. 2 for both species is , whereas since both β_s,i and t_acc,i scale with the particle species’ mass and neglecting Fermi-II processes one finds γ_max,i = (β_s,it_acc,i)^-1. Note that all these properties result from the self-consistency of our model, which is time-dependent and hence able to reproduce outbursts of blazars by replacing the steady state values of B, and/or Q_0,i by time-dependent ones. In the hard sphere approximation t_acc,i ∝ m_i ⇒ K_{| | ,i} ∝ m_i applies, i.e. the particle’s mean free path ⟨ l ⟩ ∝ m_i as expected. The dependence of t_acc on the particle’s mass is a rather simple assumption based on the gyro-frequency. A more detailed description would require in-depth knowledge of the turbulent magnetic spectrum and particle scattering properties, which is not available for AGN. Using the hard-sphere approximation as the simpified model assumption to connect the acceleration timescale to the microphysics of the jet, in principle any value of is allowed. However, we keep in mind that values t_acc,i<t_esc,i resulting in spectra harder than s = 2 are rather hard to motivate in the light of non-relativistic shock acceleration. In this case s is correlated to the compression-ratio r of the shock as s = (r + 2) / (r − 1) which may be derived from microphysical considerations incorporating bulk plasma speeds and return propabilities of the test-particle. This can be connected to the kinetic equation used (see e.g. Dermer 2013), but is not carried out as the bulk speeds of the underlying plasmoids are not feasible experimentally at this point. Hence, the power-law particle spectra generated in the acceleration zone may be identified as shock acceleration or a closely related process. We note that compression ratios r> 4 and, thus, harder particle spectra can be found, for example using non-linear simulations. In addition it should be mentioned that not using the hard-sphere approximation, i.e. t_acc,i → t_acc,i(γ_i), will yield spectra different from power-laws. This case would naturally also require a more sophisticated treatment of t_esc,i.

Fig. 2 Evolution of p+ and e− of freshly injected particles at γ0,i = 2.5 in the acceleration zone, assuming a → ∞. The separation of arises self-consistently during the simulation. p+ reaches the steady state a factor of p+/ e− slower than e−, see text for details.

As every particle escaping the acceleration zone enters the radiation zone, the particle spectra of the acceleration zone (see Fig. 2) serve as the injection function for the emission region.

Unlike in the acceleration zone, energy gains of the particles are not accounted for at this site. The particles under consideration here will cool because of synchrotron, inverse Compton, and photo-hadronic processes when they radiate. The kinetic equations for electrons, secondary positrons, and protons, as well as the radiative equation are solved self-consistently and time-dependently. Although protons are introduced merely via two new parameters in comparison to the self-Compton case described in Weidinger et al. (2010), namely the injected number density Q_0,p⁺ and the corresponding γ_0,p⁺, this has a major effect in the radiation zone, especially when time-dependency is considered, since photo-hadronic processes emerge as a significant contribution. This contribution has an impact on the radiative equation and, more importantly couples the kinetic equation of the electrons/protons to the radiative equation in a non-linear way as free pair-processes will play a major role as well, in the case of non vanishing proton densities.

For the protons the kinetic equation in the radiation zone hence yields (4)with the parameter b< 1 ensuring particle conservation. Following the discussion of Blumenthal (1970) and Berezinsky & Grigoréva (1988) the losses due to photo-meson production P_pγ can be neglected against the dominating synchrotron losses for most blazars, substantially lowering computational costs. Highly relativistic protons with γ_p⁺> 10⁵ are subject to photo-meson production with the radiation field within the emitting region, mainly pion production (5)The produced pions are unstable particles and decay into stable e^± (and γs) via the muon channel, thereby producing neutrinos of flavours (6)To calculate the production rate of stable particles from pion production we use the Kelner & Aharonian (2008) parametrization (including their erratum) of the SOPHIA Monte Carlo results (Mücke et al. 2000) for these processes. Hence, we assume the photo-hadronic interactions to be instantaneous compared to the synchrotron loss timescale, i.e. we do not account for synchrotron losses and radiation of the intermediate particles. In high magnetic fields this error remains small (Mücke et al. 2003). When considering fairly soft power-law particle spectra, as in typical blazars, Bethe-Heitler pair-production as considered in Petropoulou & Mastichiadis (2013) can be neglected since its contribution is orders of magnitude below the p⁺-synchrotron radiation, see e.g. Böttcher et al. (2013). Even though the Bethe-Heitler process sets in at lower energies, the resulting pairs in magnetic fields as low as considered here will not produce a significant contribution to the SED, although this might be different for other blazar flavours e.g. Murase et al. (2014). Nevertheless, photon quenching and runaway secondary production (Petropoulou et al. 2013) will occur in our model as well through the photo-meson channel, given the corresponding parameters (i.e. p⁺-densities).

The resulting γs from the π⁰ decay as well as synchrotron radiation of the secondary e^± from the π^± decay will partially be in the optically thick regime. Pair production (7)will occur, mainly with the synchrotron radiation field of the primary accelerated electrons. Therefore, electromagnetic cascades will emerge until the resulting synchrotron radiation is visible in the optically thin regime for pair-production. These cascades have a major non-linear effect on the kinetic equations of the electrons, protons and photons. The kinetic equations for the electrons and positrons in the radiation zone thus are No accelerated primary positrons are assumed (n_e⁺ = 0). The inverse Compton loss rate (P_IC) for isotropically distributed e^± is calculated exploiting the full Klein-Nishina cross section in Eq. (11) from e.g. Blumenthal & Gould (1970), (10)with r₀ = e²/ (mc²), hν = αmc², and q = α/ (4α₁γ²(1 − α/γ)) as the scattering parameter, i.e. α₁ being the incident photon’s energy, the scattering geometry allows for α₁<α ≤ 4α₁γ²/ (1 + 4α₁γ) which is used to determine α_max(11)The model discussed here does not assume external photon fields, the consideration of anisotropic photon fields as in Dermer et al. (2009) and Hutter & Spanier (2011) is, therefore, not necessary.

The pair-production rate in Eqs. (8) and (9) for photons of energy α₁ and α₂ respectively is calculated using the approximation (12)of Böttcher & Schlickeiser (1997) for isotropically distributed photons for the blob’s photon field N_ph with itself, with . Which has an error in energy of less than three percent over the whole electron distribution compared with the absorbed photons Eq. (18) and is numerically much more stable than the full production rate (Eq. (26) in their paper). The injection rate of secondary e^± hence is (13)with the corresponding Φ_± parameterization function of Kelner & Aharonian (2008).

The photon distribution in the radiation zone, which is eventually beamed towards the observer, reads (14)taking all mentioned processes into account. The photon escape timescale is the light crossing time. With the synchrotron emissivity R_s as the Melrose approximation (15)see Weidinger et al. (2010) and the inverse Compton production rate (16)using the full Klein-Nishina cross section Eq. (10). The photons from the π⁰ decay are calculated analogously to Eq. (13), again with the corresponding Φ₀ parametrization of the SOPHIA results: (17)In optically thick regimes the photon field is absorbed either in the low energies as a result of synchrotron self absorption by the emitting electrons/positrons themselves (α_SSA) for which the monochromatic approximation (Weidinger & Spanier 2010) is used, or in the VHEs because of e^±-pair-production. The photon annihilation coefficient for e^±-production is calculated using the exact result of Coppi & Blandford (1990) for isotropic pair plasmas, (18)with μ_max = max(−1, 1 −2x^-1), where x = αα₁ and (19)being the full cross section for γγ-pair-annihilation. The photon distribution N_ph is transformed into the observer’s frame using the beaming pattern for isotropic photon distributions in a sphere while accounting for the redshift in order to achieve the model SED of the considered blazar.

Within the model assumptions it is possible to explain the VHE peak in the typical blazar spectra either as inverse Compton upscattering of synchrotron photons or as proton synchrotron radiation consistently accompanied by cascaded radiation from the photo-hadronic interactions, just depending on the chosen parameters. We note that the p⁺-synchrotron radiation important for the occurrence of the second peak in a blazar’s SED when considering the hybrid-case naturally becomes dominant at proton densities considered for substantial photo-meson production. Thanks to the numerical treatment of the processes not relying on Monte Carlos processes, we can compute inter-band lightcurves even in the hadronic case making the potential information from outbursts of blazars accessible for physical interpretations. In the hadronic case the inter-band lightcurves exhibit more complex features (highly dependent on the chosen set of parameters) like orphan flares, typical time-lags, etc., unknown to purely leptonic models.

3. Application to 1 ES 1011+496

Despite prominent absorption lines used to pin down the redshift to z = 0.212 the blazar 1 ES 1011+496 was originally classified as a high frequency peaked BL Lac object (The MAGIC Collaboration 2007). The VHE emission of 1 ES 1011+496 was discovered with the MAGIC Air-Cherenkov telescope in 2007, and it was the most distant blazar observed in γ rays while down to the present day this is 3C 279. However, the true nature of this blazar seems to be quite unclear compared to well studied objects like Mkn 501, PKS 1218+304 or 3C 279 (The Fermi Collaboration 2011a; Weidinger & Spanier 2010; Böttcher et al. 2009). The first and only multiwavelength observation, including MAGIC in the VHEs, Swift in X-rays and KVA in the optical band was in 2008 (Reinthal et al. 2011). During this campaign the emission of 1 ES 1011+496 was more or less steady with slight spectral variability in the optical and X-rays having no counterpart in the γ rays (Reinthal et al. 2011), which of course may be due to the low observed flux and sensitivity of the instrument. The long-term optical monitoring of 1 ES 1011+496 also indicates variability consistent with synchrotron emission and magnetic field fluctuations (Böttcher et al. 2010). The multiwavelength data of 1 ES 1011+496 will now be used to study the properties of this interesting blazar in more detail with our model. The one-year Fermi-LAT butterfly is considered as an upper limit, since the MWL data represents a relatively low flux level of 1 ES 1011+496.

Fig. 3 Simultaneous data from Reinthal et al. (2011) with KVA in the optical, Swift XRT in the X-rays, and MAGIC in the VHEs as measured. Blue symbols indicate the low emission state used in the modelling. The high state inferred from the slight variability is shown in green. The grey butterfly represents the first-year catalogue spectrum of Fermi-LAT (The Fermi Collaboration 2010). The grey (online: brown) dashed curve displays the modelling attempt with a low magnetic field and no protons present in the jet (Table 1a); the black solid curve is due to the modelling with high B-field co-accelerated protons (Table 1b). The model SEDs have been EBL-absorbed using the model of Primack et al. (2005).

3.1. The spectral energy distribution

In Fig. 3 the measured multiwavelength spectrum from Reinthal et al. (2011) and two model SEDs are shown. We note that we use the EBL absorbed data points in this paper. The parameters of the models can be found in Tables 1a and b. The rather hard spectrum in the optical combined with a photon index of s = 2.32 in the Swift XRT band makes it difficult to explain the SED of 1 ES 1011+496 in terms of a self-consistent SSC-Ansatz already in the steady state (brown dashed line in Fig. 3). The best parameters in the self-Compton limit, i.e. a low magnetic field unable to confine non-thermal protons within the considered region, are shown in Table 1a. A spectral index of 2.2 for the electrons with a magnetic field of 0.18 G can neither reproduce the narrow synchrotron peak nor the spectral properties when computed self-consistently. Of course one can understand the broadband SED of 1 ES 1011+496 with a purely leptonic model, see e.g. Reinthal et al. (2011) itself, but in those cases the underlying electron spectra do not arise in the modelling process but are adhoc assumptions which spectral indices and breaks requiring physical motivation put in by hand. Consistent cooling breaks fail to reproduce the synchrotron peak at low B-fields, automatically yielding a relatively high magnetic field being present within the jet. In such a case r_gyr ≪ R_rad is fulfilled even for relativistic protons making them effective emitters within the jet, hence being co-accelerated in the acceleration zone, leading to a hybrid spectrum of this particular blazar, see Fig. 3. The parameters used can be found in Table 1b. In a magnetic field of 8.0 G the power-law electrons (with the corresponding t_acc) cannot exceed Lorentzfactors of γ_e ≈ 2.5 × 10⁴ because of synchrotron losses, which together with the relatively high γ_0,e⁻ leads to a narrow peaked νF_ν-spectrum in the X-ray band, explaining the SED of 1 ES 1011+496 there. The Fermi-LAT data on other blazars also suggests that a high γ_0e⁻ is present, e.g. Reinthal et al. (2011); Weidinger et al. (2010). These values of γ_0i are expected since Fermi acceleration requires pre-accelerated, non-thermal particles to be efficient. Hence the ratio of γ_0p⁺/γ_0e⁻ depends on this process and might give hints on the nature of this mechanism, e.g. trapping of low energetic e⁻. In magnetic fields as high as inferred here, inverse Compton scattering of synchrotron photons is negligible. The second peak in the spectrum thus consists of synchrotron radiation of highly relativistic protons, consistently accelerated to γ≈ 8.7 × 10¹⁰ in the model, as well as pair cascades arising from photo-hadronic interactions of these protons with the synchrotron photons. With the parameters shown in Table 1b the measured MWL spectrum of 1 ES 1011+469 is explained well, introducing no new parameters other than the number of injected low energy protons and their injection energy (or the electron/proton ratio and – energy ratio) as compared to the (simple) SSC-limit. The injected luminosity (i.e. the jet luminosity as the blob moves along) is L_p = 3.0 × 10⁴² erg s^-1 in hadrons, assuming a vanishing angle to the line of sight. This is significantly below the Eddington limit for the black hole assumed to be present in 1 ES 1011+496 (Fan et al. 2007). Proton spectra extending down to γ_0p⁺ = 1 would dramatically increase the kinetic luminosity without significant radiative output. The photon indices of the second peak are not mapped to the synchrotron spectrum of the radiating protons directly, since it also contains redistributed radiation (see Fig. 4) hardening the VHE peak. The spectral indices of the underlying particle species (s = 2.25 for both) is consistent with diffusive shock acceleration at a strong shock.

Fig. 4 Intrinsic SEDs (unaffected by EBL absorption) of 1 ES 1011+496 (Fig. 3) as inferred by the model, parameters in Table 1. The solid black line shows the hybrid SED with the individual components: proton synchrotron emission (dashed black line) and redistributed synchrotron emission by e±-pairs initiated by photo-hadronic interactions (dash-dotted black line) only observable through pair-cascades. The direct contribution of stable end-products from photo-pion production are negligible in the observable regime peaking above 1027 Hz opaque to γγ-pair-production, see solid brown line (π+ → e+; other processes are below 10-13 erg cm-2 s-1) and only affect the SED indirectly via γγ-pair-production of the nth generation.

3.2. Intrinsic spectrum of 1 ES 1011+496

To resolve the second peak in the spectrum of 1 ES 1011+496 it is much more convenient to take a look at the intrinsic SED unaffected by the EBL, shown in Fig. 4. As already mentioned the first peak in the SED consists of synchrotron radiation of the primary electrons within the radiation zone. The situation for the second hump, however, is slightly more complex. The relevant contributions to the second peak of 1 ES 1011+496 are proton synchrotron photons of the highly relativistic primary p⁺ with Lorentzfactors up to γ_p ≈ 10¹⁰ and cascade radiation. The latter is synchrotron radiation of the stable products arising from photo-meson production which is redistributed from the optically thick regime for pair-production until observable. As one can infer from Fig. 4 (solid brown line) the direct contribution of pγ interactions is negligible thanks to the dominance of the proton synchrotron peak. That is also why Bethe-Heitler pair-production with a lower threshold than photo-hadronic processes can be neglected against the proton synchrotron emission in magnetic fields of G) required to confine the protons within the emitting region of a typical blazar. From Fig. 4 it is also clear that the maxima of the first generation radiation of e^± and γs from pγ-interactions are above 10²⁸ Hz.

It should be noted that the attenuation of the spectrum by the EBL may have secondary effects. It has been discussed in the literature (cf. Neronov & Vovk 2010) that EBL absorption yields electron-positron pairs, which can upscatter CMB photons to GeV energies. The detection or non-detection of these GeV photons by Fermi, can give rise to the determination of the intergalactic magnetic field (IGMF; assuming that diffusion of electron-positron pairs in the IGMF is the only process preventing the Compton scattering). The integrated flux above 1 TeV in the hybrid scenario is an order of magnitude higher than the hard spectra assumed in Neronov & Vovk (2010) and as the redshift of 1 ES 1011+496 is comparable, we find a crude lower bound for the IGMF of G), which has been reported in other analyses, e.g. Tavecchio et al. (2010). This is an extremely rough estimate, since different regions of the EBL are touched in the hybrid case compared to purely leptonic spectra which cut off at lower energies, additionally leading to deviating IC cooling timescales. It should not be neglected here that plasma effects (namely the pair-beam instability) with substentially smaller mean free paths are passionatly discussed as the dominating cooling process for TeV beams, not producing a radiative signature in the Fermi-LAT band, even in the absence of an IGMF (Schlickeiser et al. 2013; Broderick et al. 2012). Furthermore the pair creation probability inside the jet, but outside the actual source region has not been modelled, which may substantially alter the TeV spectrum of a blazar injected into the intergalactic medium.

3.3. Variability

Yet another advantage of the numerical approach is that it is not only possible to use exact cross sections including all non-linear interactions, but that a time-dependent treatment is feasible and variability becomes accessible even in the lepto-hadronic case. The structure of the VHE peak in the hybrid scenario including leptonic and hadronic emission will involve various timescales of variability which can be investigated with our model. In this section we introduce a possible flaring scenario for 1 ES 1011+496, where more primary e⁻ and p⁺ are injected into the steady state acceleration zone for a certain amount of time Δt, i.e. (20)This scenario might occur within the jet of a blazar as the blob moves through strong density fluctuations along the axis of the outflow. Unlike the SSC models, we expect long time-lags and different timescales to occur in the multiband lightcurve serving as a hadronic fingerprint of the considered blazar. All timescales and times are given in the observer’s frame.

Figure 5 shows the response in the optical (blue), X-rays (red), and gamma rays (black) when (arbitrarily) setting Δt = 3.86 h and x_e = 2.5, x_p = 8.5. Of course these parameters need to be interpreted and will be different observing an actual outburst of 1 ES 1011+496 in various energy bands. The definition of the energy ranges in the multiband lightcurves of Fig. 5 can be found in Fig. 6. A flare occurs almost instantaneously in the optical (2 × 10¹⁴ Hz to 9 × 10¹⁴ Hz) and X-ray (1 × 10¹⁷ Hz to 9 × 10¹⁷ Hz) energy ranges due to the higher primary electron density accelerated to high energies. Typical timescales of this first outburst are the acceleration and synchrotron loss timescale of electrons, see inset in Fig. 5. Along with the flare in the optical comes an outburst at VHE (1 × 10²⁵ Hz to 2 × 10²⁶ Hz) with a slightly different rising behaviour due to enhanced reprocessed radiation of the initial photo-hadronic interactions, since more target (synchrotron) photons are provided during the outburst. While the electron synchrotron peak shows spectral behaviour similar to SSC modelling, the photon index in the VHE softens as the cascaded radiation dominates over the proton synchrotron emission during the first flare, see Fig. 6. After approximately 0.7 days the emission initiated by the enhanced Q_0,e⁻ has cooled down to its steady-state value. The additionally injected protons (alongside the electrons) need to be accelerated with t_acc,p to high γs before emitting synchrotron radiation efficiently. This causes the delay t_delay between the first, MWL flare and the second, orphan flare in the VHE due to p⁺-synchrotron radiation with its typical timescales, yet another feature not expected in a purely leptonic model. In contrast to the first outburst the photon index in the VHE hardens (inset of Fig. 6). Both are initiated by one single enhanced injection into the numerically considered region. Hence the VHE exhibit both characteristic timescales, those of electrons and those of protons in the magnetic field of the blob, and it is insufficient to limit theoretical investigation to p⁺-timescales only.

Fig. 5 Multiband lightcurve during the assumed scenario of density fluctuations along the jet axis explained in the text. The inset shows the definition of the snapshots during the first outburst in Fig. 6. The timescale of this flare corresponds to the synchrotron timescale of the electrons with a slight delay in the gamma rays caused by the cascades that need to be built up first. With a delay of about 4.28 days a second, orphan flare occurs in the VHEs when the freshly injected protons have been accelerated to high Lorentz factors to radiate efficiently. The timescale here corresponds to the synchrotron loss timescale of protons.

Fig. 6 Spectral evolution (18.5 min snapshots) during the imposed flaring scenario of 1 ES 1011+496 (intrinsic spectrum). The VHE peak rises with the synchrotron peak because of the enhanced seed photons for hadronic interactions and the redistribution via γγ → e±. The inset shows the second, orphan flare in the gamma rays (see also Fig. 5) due to enhanced proton synchrotron radiation occurring with a delay, since the protons need time to accelerate to high energies.

Fig. 7 Spectral evolution (43.2 min snapshots) during the imposed flaring scenario by setting xe = 2.5 in the leptonic case of 1 ES 1011+496 (intrinsic spectrum). The SSC limit shows the simple flaring behaviour as expected, without spectral evolution in the VHE and orphan flares do not occur. The inset shows the corresponding lightcurves in the X-rays and gamma rays exhibiting a slight delay caused by the IC upscattering of the enhanced synchrotron photons.

In Fig. 7 a hypothetical outburst with x_e = 2.5 and Δt = 3.86 h in terms of the purely leptonic model (Table 1a) is displayed for comparison, the inset shows the corresponding lightcurves. The behaviour is quite straightforward with its possible soft- and hard lags occurring, in contrast to the relatively complex time-dependent behaviour in the hybrid model. There is no noteworthy spectral evolution in the VHE and orphan flares do not occur. The timescales as shown in the lightcurves mainly differ from those of Table 1b) as a result of the significantly lower magnetic field and are consistent with the modelling of comparable blazars (Weidinger & Spanier 2010).

4. Discussion

We have presented and applied a fully self-consistent and time-dependent hybrid emission model for blazar jets, taking all relevant processes into account, namely acceleration and synchrotron emission of electrons and protons, inverse Compton scattering, and intrinsic photo-hadronic interactions as well as γγ-pair production. The model includes the self-consistent SSC limit by setting Q_0,p⁺ → 0, see e.g. Weidinger et al. (2010) for a detailed description. However magnetic fields of G) will confine relativistic protons within the emitting region with a typical size of cm) to cm) as inferred from variability and causality (Weidinger et al. 2010; Böttcher & Dermer 1998). The numerical approach allows for a detailed treatment of the individual processes (see Chap. 2), including all non-linearities arising from the coupling of the photons to the leptons via photo-meson production and the creation of e^±-pairs with different timescales. Hence the model spectra and lightcurves in case of outbursts become far more complex than in SSC models, even if non-linear cooling is included (Weidinger & Spanier 2010; Zacharias & Schlickeiser 2013). Not only do orphan flares occur, but also variability associated with the different timescales for p⁺ and e^± in the system. More complex behaviour will arise if photon-quenching becomes effective, in the cases when the electromagnetic cascades serve as targets for photo-meson production. It is obvious that photo-hadronic cooling effects must not be neglected in those cases, but even without runaway production there is vast, qualitatively different, variable multiband behaviour, delicately depending on the chosen set of parameters, not only the behaviour 1 ES 1011+496 exhibits in our example. In addition typical time-lags are directly bound to the timescales, e.g. the acceleration timescale, in the system. Hence these appealing features could be used to determine typical acceleration rates in astrophysical jets, otherwise inaccessible. With future experiments a systematic approach of multiband variable behaviour is at hand, leading towards an investigation of the energy dependence of the involved timescales comparing lags and durations in different energy bands in our model, leading to advances in e.g. relativistic acceleration of particles, still not understood completely (Bykov et al. 2012; Ellison et al. 2013; Sironi et al. 2013).

1 ES 1011+496 was used as an example blazar with its low state emission already pointing out a relatively high magnetic field to be present within its jet, see Chap. 3. The steady state observed in the MWL campaign including the MAGIC telescopes and the Swift satellite is reliably described with our hybrid model, while a consistent leptonic approach fails to model the synchrotron peak of 1 ES 1011+496 (see Table 1 for the parameters used). One can of course assume various breaks with corresponding spectral indices in the electron distribution to describe the X-ray peak leptonically, but in general a physical motivation for these is rather hard to find. Hence, self-consistently speaking, the VHE peak is not due to inverse Compton photons, but consists of proton synchrotron radiation and emerging cascade radiation. We note that merely two additional parameters need to be introduced in the complex lepto-hadronic model with respect to the Weidinger et al. (2010) or Tavecchio et al. (2001) SSC models. The properties of the underlying electron and proton densities, like the gap between γ_max,e and γ_max,p or spectral indices and breaks, arise consistently as a result of acceleration and cooling during the modelling process.

This has a certain impact on the variability observable from 1 ES 1011+496, investigated using an imposed flaring scenario simulating density fluctuations along the jet-axis (Eq. (20)). Quite obviously the timescales connected with e⁻-synchrotron emission is lower than in the typical SSC-case where the magnetic fields are low. More importantly, one can observe patterns in the multiband lightcurves, see Fig. 5. In the imposed scenario, the first flare of 1 ES 1011+496 covers all considered energy bands with the timescale of the electrons, even in the VHEs, since it is caused by the freshly injected e⁻ and the enhanced cascade radiation. With a time-lag of approximately 4.28 days a second, orphan flare occurs in the VHEs because of p⁺-synchrotron radiation with its typical timescale. It took the protons (with their higher t_acc,p) this time to be accelerated to the energies where they can emit efficiently. The fact that leptonic timescales are involved in the VHE emission as well, allows short-time variability in this region even in the hybrid case. This is not necessarily true in every hadronically dominated blazar, see e.g. Böttcher et al. (2009) where p⁺-synchrotron radiation dominates over the cascades. In comparison with the outburst in the purely leptonic case, one can see that strong spectral features in the VHE during the flare are only present within the hybrid model with significant hardening and softening.

The multiband lightcurves can therefore be used as a hadronic fingerprint of an individual blazar by identifying the characteristic patterns using the presented model. The tangible patterns will be different in each blazar according to the physical parameters present within the jets, see e.g. Weidinger & Spanier (2011), but typical timescales and timelags can only be found in the hybrid case and thus be used to identify blazars accelerating both electrons and hadrons without relying on the detection of neutrinos of an individual source, e.g. with IceCube, which is quite difficult (Abbasi et al. 2011; Mücke & Protheroe 2001; Vissani et al. 2011; Vitells & Gross 2011). Consequently this paper suggests more long-term MWL observations of individual sources as done with Fermi-LAT and VERITAS of Mkn 501 or Mkn 421 (The Fermi Collaboration 2011a,b) to reveal those responsible for the recently detected extragalactic high energy neutrinos (Aartsen et al. 2013). Of course, AGN accelerating protons would be good candidates of ultra high energy cosmic ray production sites among blazars, since proton energies up to E_p = 10²¹ eV are reached in the observer’s frame, and it is thus important to identify the subset of those actually confining non-thermal protons, a task which can be addressed by our model. Systematic modelling of all types of blazars with one single model could also help to understand the differences and commonalities of each blazar flavour (FSQR, IBL, LBL, and HBL), because the parameters are obtained during the modelling. They are not biased using a specific Ansatz like leptonic or hadronic in the first place when producing the model spectrum. Thus, models like ours are crucial tools that can be used to interpret data in times of extensive multi-messenger campaigns and stacking astrophysics with the upcoming CTA, IceCube detecting first PeV neutrinos, and cosmic ray detectors like AUGER and HAWC in the future.

Acknowledgments

We thank the referee for a helpful and detailed report. M.W. wants to thank Elitenetzwerk Bayern and the Graduate School 1147 for their support, F.S. acknowledged support from National Research Foundation through a “Multiwavelength” grant

References

All Tables

All Figures

	Fig. 1 Basic model geometry. All escaping particles of the acceleration zone serve as injection for the radiation zone, i = e± for electrons/positrons and i = p+ for protons.
In the text

	Fig. 2 Evolution of p+ and e− of freshly injected particles at γ0,i = 2.5 in the acceleration zone, assuming a → ∞. The separation of arises self-consistently during the simulation. p+ reaches the steady state a factor of p+/ e− slower than e−, see text for details.
In the text

Fig. 3 Simultaneous data from Reinthal et al. (2011) with KVA in the optical, Swift XRT in the X-rays, and MAGIC in the VHEs as measured. Blue symbols indicate the low emission state used in the modelling. The high state inferred from the slight variability is shown in green. The grey butterfly represents the first-year catalogue spectrum of Fermi-LAT (The Fermi Collaboration 2010). The grey (online: brown) dashed curve displays the modelling attempt with a low magnetic field and no protons present in the jet (Table 1a); the black solid curve is due to the modelling with high B-field co-accelerated protons (Table 1b). The model SEDs have been EBL-absorbed using the model of Primack et al. (2005).

In the text

Fig. 4 Intrinsic SEDs (unaffected by EBL absorption) of 1 ES 1011+496 (Fig. 3) as inferred by the model, parameters in Table 1. The solid black line shows the hybrid SED with the individual components: proton synchrotron emission (dashed black line) and redistributed synchrotron emission by e±-pairs initiated by photo-hadronic interactions (dash-dotted black line) only observable through pair-cascades. The direct contribution of stable end-products from photo-pion production are negligible in the observable regime peaking above 1027 Hz opaque to γγ-pair-production, see solid brown line (π+ → e+; other processes are below 10-13 erg cm-2 s-1) and only affect the SED indirectly via γγ-pair-production of the nth generation.

In the text

Fig. 5 Multiband lightcurve during the assumed scenario of density fluctuations along the jet axis explained in the text. The inset shows the definition of the snapshots during the first outburst in Fig. 6. The timescale of this flare corresponds to the synchrotron timescale of the electrons with a slight delay in the gamma rays caused by the cascades that need to be built up first. With a delay of about 4.28 days a second, orphan flare occurs in the VHEs when the freshly injected protons have been accelerated to high Lorentz factors to radiate efficiently. The timescale here corresponds to the synchrotron loss timescale of protons.

In the text

Fig. 6 Spectral evolution (18.5 min snapshots) during the imposed flaring scenario of 1 ES 1011+496 (intrinsic spectrum). The VHE peak rises with the synchrotron peak because of the enhanced seed photons for hadronic interactions and the redistribution via γγ → e±. The inset shows the second, orphan flare in the gamma rays (see also Fig. 5) due to enhanced proton synchrotron radiation occurring with a delay, since the protons need time to accelerate to high energies.

In the text

Fig. 7 Spectral evolution (43.2 min snapshots) during the imposed flaring scenario by setting xe = 2.5 in the leptonic case of 1 ES 1011+496 (intrinsic spectrum). The SSC limit shows the simple flaring behaviour as expected, without spectral evolution in the VHE and orphan flares do not occur. The inset shows the corresponding lightcurves in the X-rays and gamma rays exhibiting a slight delay caused by the IC upscattering of the enhanced synchrotron photons.

In the text