Hadronic supercriticality as a trigger for GRB emission
Abstract
We explore a onezone hadronic model that may be able to reproduce ray burst (GRB) prompt emission with a minimum of free parameters. Assuming only that GRBs are efficient highenergy proton accelerators and without the presence of an ab initio photon field, we investigate the conditions under which the system becomes supercritical, i.e. there is a fast, nonlinear transfer of energy from protons to secondary particles initiated by the spontaneous quenching of protonproduced rays. We first show analytically that the transition to supercriticality occurs whenever the proton injection compactness exceeds a critical value, which favours high proton injection luminosities and a wide range of bulk Lorentz factors. The properties of supercriticality are then studied with a timedependent numerical code that solves concurrently the coupled equations of proton, photon, electron, neutron and neutrino distributions. For conditions that drive the system deep into the supercriticality we find that the photon spectra obtain a Bandlike shape due to Comptonization by cooled pairs and that the energy transfer efficiency from protons to rays and neutrinos is high reaching . Although some questions concerning its full adaptability to the GRB prompt emission remain open, supercriticality is found to be a promising process in that regard.
keywords:
astroparticle physics – instabilities – radiation mechanisms: nonthermal –gamma ray burst: general1 Introduction
The prompt emission of gammaray bursts (GRBs) is observed in the keV MeV energy band (Preece et al., 2000; von Kienlin et al., 2014), with fluences that generally range between erg cm and erg cm, where the lower limit does not necessarily reflect an intrinsic property of GRBs but depends on the sensitivity of the detectors. The fluence distribution of bursts detected by the Fermi/GBM peaks at erg cm (von Kienlin et al., 2014), which is considered as the typical GRB fluence. GRB light curves are highly variable and have a complex structure: they consist of several pulses (10100) each of them having typical width mss (Norris et al., 1996; Nakar & Piran, 2002), thus making the total duration of the burst longer, e.g. s. The GRB spectra can, in most cases, be described by a smoothly connected broken powerlaw (Band et al., 1993, 2009) with values of the break energy clustered around MeV in the observer’s frame, while even higher peak energies ( MeV) have been detected (e.g. Goldstein et al. 2012). The typical photon indices of the spectrum below and above the peak are and , respectively (Preece et al., 2000; Goldstein et al., 2012).
The origin of the GRB emission is still an open issue, although various models have been proposed during the past decades trying to address all or most of the above properties. Borrowing the terminology from the field of blazar modelling (see e.g. Böttcher 2007; Boettcher 2010, for reviews) GRB emission models can be divided in two categories, namely leptonic and hadronic, according to the species of the radiating particles.
The former try to attribute the gammaray emission by employing radiation processes of relativistic electrons (and/or positrons). The classical scenario of the prompt emission, which belongs to the first category, is the optically thin synchrotron model (e.g. Katz 1994; Sari et al. 1996; Tavani 1996; Chiang & Dermer 1999), where the kinetic energy of the flow is dissipated via shocks and the prompt emission is the result of synchrotron radiation of relativistic electrons. The difficulties that this scenario has in dealing with several issues, such as the lowenergy photon index (Crider et al., 1997; Preece et al., 1998, 2000) and the fine tuning of parameters required to explain the MeV peak of GRB spectra (see also Beloborodov 2010 for relevant discussion), motivated works on alternative scenarios. For example, variants of nonthermal emission models were discussed in order to overcome the ‘lineofdeath’ problem: (i) effects of adiabatic and/or inverse Compton (IC) cooling in the KleinNishina regime on the low energy part of the synchrotron spectrum (Derishev et al., 2001; Wang et al., 2009; Daigne et al., 2011); (ii) synchrotron emission from an electron distribution with a smooth low energy cutoff and an anistropic distribution in pitch angles (Lloyd & Petrosian, 2000); (iii) the emission observed in the BATSE energy band 20 keV1 MeV was interpeted as the result of IC scattering of slow cooling electrons on the selfabsorbed part of the synchrotron spectrum (Panaitescu & Mészáros, 2000); (iv) synchrotron selfCompton emission under the assumption of continuous electron acceleration (Stern & Poutanen, 2004); (v) jitter radiation emitted by relativistic electrons moving in non uniform small scale magnetic fields (Medvedev, 2000); (vi) synchrotron electron cooling in a decaying magnetic field (Pe’er & Zhang, 2006); (vii) gammaray emission through the Comptondrag process (Lazzati et al., 2000). The socalled photospheric models, where the radiation is released when the outflow becomes transparent, constitute an interesting alternative to the nonthermal ones. If the energy is dissipated at the very inner parts of the outflow, it thermalizes and the radiation that escapes from the GRB photosphere has a quasithermal spectrum that peaks at MeV (e.g. Goodman 1986; Thompson 1994; Beloborodov 2010). In the presence of continuous energy dissipation, however, the resulting spectra may obtain a nonthermal appearance via Comptonization of the quasithermal emission by thermal electrons (Mészáros & Rees, 2000; Pe’er et al., 2006; Giannios, 2006, 2012).
Hadronic models for the GRB prompt emission constitute a viable alternative. They are built upon the common basis that the (sub)MeV ray emission, which serves as the target field for photopion interactions, is not of hadronic origin but it is either the synchrotron radiation of primary electrons (e.g. Dermer & Atoyan 2003; Asano & Inoue 2007; Murase 2008) or the emission from the photosphere itself (e.g. Gao et al. 2012; Asano & Mészáros 2013). The highenergy part of the gammaray spectrum (100 MeV) is typically explained by relativistic proton synchrotron radiation (Vietri, 1997; Totani, 1998) or by protoninduced cascades (Dermer & Atoyan, 2006; Asano & Inoue, 2007; Asano et al., 2009). The latter scenario has been applied to explain the underlying powerlaw components seen in some bright Fermi bursts (e.g. GRB 090902B (Abdo et al., 2009); GRB 080319B (Racusin et al., 2008)), which extend from the hard Xrays up to GeV energies and do not agree with simple extrapolations of the MeV spectrum (Asano et al., 2010). In any case, the suggestion that GRBs are the sources of ultrahigh energy cosmic rays (UHECRs) (Waxman, 1995b; Vietri, 1995; Murase et al., 2008) makes hadronic models attractive. Moreover, the associated highenergy neutrino emission has been calculated in various studies (Waxman & Bahcall, 1997; Murase, 2008; Gao et al., 2012; He et al., 2012; Zhang & Kumar, 2013; Asano & Meszaros, 2014; Baerwald et al., 2014; Reynoso, 2014; Petropoulou et al., 2014) and now starts becoming testable by ongoing observations (IceCube Collaboration, 2013; Aartsen et al., 2014b, a).
However, onezone hadronic models are inherently more complex than pure leptonic ones, since they require modelling of the coupled emission and energy loss processes (see e.g. Dimitrakoudis et al. 2012a) between various species in order to track the evolution of the different components (protons, neutrons, pairs, mesons, neutrinos, photons). Several of these feedback processes were proven to give rise to radiative instabilities (e.g. Stern & Svensson 1991; Kirk & Mastichiadis 1992; Mastichiadis et al. 2005) that share a common feature: the abrupt, i.e. in a few dynamical times, release of energy that is initially stored in protons and is subsequently transfered to photons. The proton synchrotron pairproduction instability for example was proposed to give rise to gammaray emission peaking at MeV (e.g. Kazanas et al. 2002; Mastichiadis & Kazanas 2006) offering, at the same time, a physical connection between the prompt and afterglow phases (Mastichiadis & Kazanas, 2009; Sultana et al., 2013). Moreover, possible implications of the automatic ray quenching instability (Stawarz & Kirk, 2007; Petropoulou & Mastichiadis, 2011) were studied in Petropoulou & Mastichiadis (2012a, b) in the context of hadronic blazar emission.
In the present work we extend this analysis by exploring its emission signatures in gammarays, highenergy neutrinos and cosmicrays for parameters relevant to GRB sources. We begin with the sole assumption of a source that is an efficient ultra highenergy (UHE) proton accelerator and sufficiently magnetized in order to confine the accelerated protons. Contrary to the majority of GRB hadronic models we do not assume an external source of photons. We show that for low values of the proton injection luminosity, the only photon field present in the source is the one emitted by the proton component mainly through synchrotron radiation. In this case, the emitted spectrum cannot be assigned to that of a typical GRB because of its spectral shape and its low luminosity. We follow the evolution of protons by balancing their losses to the respective gains of their secondaries as done in Mastichiadis & Kirk (1995) and show that if their injection luminosity exceeds a critical value, the system undergoes a transition that is triggered by the instability of automatic ray quenching. The transition is easily identified by an abrupt increase of the photon luminosity that causes the source to enter in a high photon compactness state making the energy exchange between leptons and photons dominant. This results in photon spectral shapes that in general resemble GRBs in the sense that (i) they match the required luminosity and (ii) they can be fit by a Band function. It is this selfconsistently produced radiation field that becomes the target for photopion interactions and efficiently drains energy from UHE protons, part of which is transfered to highenergy electron and muon neutrinos produced through the charged pion decay. This constitutes one of the fundamental differences between the present study and others where an ad hoc Band spectrum is assumed (see also Murase et al. 2008 for a similar approach to ours). Finally, relativistic neutrons produced from the same photopion interactions provide an effective means for UHECR escape from the source, since they are not magnetically confined and their decay time is long enough as to allow them to escape freely before converting into protons through decay. We show that for modest values of the bulk Lorentz factor () all three components, namely neutrinos, UHECRs and photons, are energetically similar.
The present work is structured as follows. In §2 we present the physical conditions of the model and derive an analytical expression for the critical proton injection luminosity. In §3.1 we present the numerical code and in §3.2 we continue with a presentation of the photon, UHECR and neutrino emission spectra; we also discuss the effects of a lower value for the highenergy cutoff of the proton distribution. We discuss our results in §4 and conclude with a summary in §5. Throughout this study we use km Mpc s, , and as an indicative value for the redshift of the gammaray source. We also introduce the notation .
2 Analytical approach
2.1 Physical conditions
The GRB lightcurve in the soft ray band (10 MeV) is variable and consists of several pulses with durations in the range s (Norris et al., 1996; Nakar & Piran, 2002). Here we adopt the internal shock scenario (see Piran (2004); Zhang & Mészáros (2004), for reviews) according to which the emitting region that corresponds to each individual pulse is modelled as a homogeneous shell with Lorentz factor that forms at a distance from the central engine. In the comoving frame the width of the shell is
(1) 
As long as the beaming angle is smaller than the opening angle of the jet, which holds during the internal shock phase, we may treat the emission region as a spherical blob of radius .
Since the typical isotropic energy emitted in rays is erg (Bloom et al., 2003; Kocevski & Butler, 2008), the isotropic ray luminosity defined as ranges between and erg/s, for a fiducial burst duration s. A minimal requirement is that , where where is the total power of the jet and equals to with , and being the kinetic, Poynting and proton^{1}^{1}1 refers to the luminosity of a powerlaw proton distribution. luminosities, respectively. Although electron acceleration at high energies is expected to take place too, here, in our attempt to keep the number of free parameters as low as possible, we assume that the injection luminosity of primary relativistic electrons is much lower than that of protons, making their contribution to the energetics and the overall spectra negligible.
We introduce next the parameters and . We use throughout the text and as indicative values, unless stated otherwise. Using the definition of the Poynting luminosity
(2) 
we may write the magnetic field strength measured in the comoving frame as follows:
(3) 
As we show next, all physical quantities in addition to and may be expressed through five essential variables: , and .
We proceed with the derivation of the respective expressions for the two basic quantities that describe the proton distribution, namely its injection compactness and its highenergy cutoff. The former, is defined as
(4) 
and using eq. (1) it is also written as
(5) 
We assume that protons are being injected into the blob after having been accelerated into a powerlaw distribution with index starting from up to a Lorentz factor which is usually determined by the balance between the acceleration and the energy loss processes. Because protons do not, in principal, suffer as severe radiative losses as electrons do, the acceleration of protons to large in GRBs is possible (e.g. Asano et al. 2009; Murase et al. 2012). For the lower cutoff of the proton distribution and the powerlaw index we use the indicative values and respectively. We note that the exact value of does not alter the main results of this work, while steeper proton spectra would increase the energy demands; for this reason we will not consider such cases here. Assuming that the synchrotron losses are the dominant energy loss process for high energy protons (see also Petropoulou et al. 2014) and that the acceleration process operates close to the Bohm diffusion limit (e.g. Giannios 2010), the typical energy loss and acceleration timescales are given respectively by and , where and . Using the above expressions and eq. (3) we find
(6) 
Thus, protons can in principle be accelerated to ultrahigh energies (UHE), i.e. eV in the comoving frame. If the gyroradius of these highly energetic protons is, however, larger than the typical size of the emission region , they cannot be confined and escape from it. Thus, the maximum energy of protons in the emission region is given by , where is derived using the Hillas criterion (Hillas, 1984), namely . Using eqs. (1) and (3) this is also written as
(7) 
Combining eqs. (6) and (7) we find that , unless
(8) 
For the fiducial parameter values used here and for the purposes of the analytical treatment presented in §2.2. it is sufficient to define the maximum proton Lorentz factor through eq. (7). However, in §3 where we study the problem numerically for different parameter sets we use the appropriate expression for .
2.2 Transition to supercriticality
The injected protons will emit synchrotron radiation that peaks in the comoving frame at , where , G, and ; here, is the typical Lorentz factor of protons that cool due to synchrotron losses within the dynamical timescale and it is given by
(9) 
Using eqs. (3) and (7) we find that
(10) 
where for simplicity we assumed (for the validity of the assumption, see Appendix A). Thus, for the fiducial parameter values used here the peak of the proton synchrotron spectrum falls well within the ray energy band. As long as
(11)  
(12) 
is satisfied^{2}^{2}2For the derivation we used the feedback criterion of automatic photon quenching – see Petropoulou & Mastichiadis (2011) and eq. (2) therein., photons at the highenergy part of the proton synchrotron spectrum can be successible to the instability of spontaneous gammaray quenching (Stawarz & Kirk, 2007; Petropoulou & Mastichiadis, 2011). According to this, the ray compactness () produced in a highly magnetized region cannot become arbitrarily high. Whenever it exceeds a critical value () it is spontaneously absorbed producing relativistic pairs that cool by emitting a large^{3}^{3}3An electron with energy will emit photons, where . number of synchrotron photons, thus providing more targets for further absorption. As the ray compactess in our framework is related to the proton injection compactess, the existence of an upper limit to is also translated to a limiting value for the proton injection compactess (). This has been already pointed out by Petropoulou & Mastichiadis (2012a), although in a different context.
We refer the reader to Appendix A for the derivation of and here we present the final result
(13) 
where the first branch is relevant for and the second otherwise. From this point on we will refer to cases with as ‘subcritical’ and ‘supercritical’ otherwise. We will also drop the logarithmic dependance and set instead the numerical factor in the bracket of eq. (13) to 20, i.e. .
Note that the above expression is valid as long as the gammarays that are spontaneously absorbed and initiate the instability are the result of proton synchrotron radiation. If the parameters are such as to push the peak of the proton synchrotron radiation to GeV energies, e.g. , then the emission from pairs produced by BetheHeitler and/or photopion interactions of protons with their own synchrotron radiation dominates in TeV energies. It can be shown that even in such cases the instability can still operate (Petropoulou & Mastichiadis, 2012a). Although the derivation of an expression similar to eq. (13) in this case is out of the scope of the present study, we will present a detailed numerical example in §3.
The absolute value of is not so important by itself. More important for the evolution of the system is the ratio which measures how deep in the supercritical regime the system is driven for given physical conditions. Using eqs. (5) and (13) this is written
(14) 
where we used the first branch of eq. (13) for simplifying reasons. In any case, both expressions of are similar. Interestingly, eq. (14) shows that the condition is satisfied for a wide range of parameter values, with lower values of being prefered. Because of the strong dependance of the ratio on , slightly different values of the bulk Lorentz factor lead to very different photon and neutrino spectra as we show in §3. It is useful, therefore, to define a ‘critical’ value of the Lorentz factor too. Setting we find
(15) 
which has very weak dependance on the parameters.
The above are summarized in Fig. 1 where we plot (solid lines) and (dashed lines) for three values of marked on the plot. Other parameters used are: , and s. The dashed lines are plotted up to (see eq. (12)), as the derived expression for the critical compactness is not relevant for larger values of . For each value of we plotted with thick lines the part of the curve that lies above . Finally, the red lines denote the loci of points with and 10. The parameter space that leads to supercriticality becomes wider as increases. For high enough values, e.g. erg/s, supercriticality is ensured for almost all values of relevant to GRBs. Note also that low values of the bulk Lorentz factor not only favour the transition to supercriticality but also correspond to .
The condition implies that the system lies deep in the supercritical regime and, as we will show in the next section with detailed numerical examples, this has the following implications: (i) the neutrino production efficiency is high, (ii) UHE protons cool down effectively through photopair and photopion processes, and (iii) the photon spectrum may be adequately described by a Band function (Band et al., 1993, 2009).
2.3 Comparison with the WB model
Waxman and Bahcall derived an elegant expression for the energy lost by protons through pion production within a dynamical timescale (see eq. (4) in Waxman & Bahcall 1997), which simply depends on the observed peak energy and ray luminosity as well as on and . Since proton cooling due to photopion interactions becomes more inefficient for larger values of , their expression can also be seen as an upper limit for
(16) 
In other words, in the WB model, pion and neutrino production is efficient for , which requires either low values of the Lorentz factor if s or modest values of , e.g. , if the ray variability is extremely fast, i.e. ms (see e.g. Waxman & Bahcall 1997; Guetta et al. 2004; Abbasi et al. 2010). In our framework, however, efficient pion production is ensured, even without the requirement of an ab initio target photon field, for parameters leading to the supercritical regime (see e.g. Fig. 9 in §3). In the previous paragraph we showed by analytical menas that the transition to supercriticality occurs for , where is defined in eq. (15). A comparison between the two Lorentz factors is shown in Fig. 2, where and are plotted with dashed and solid lines, respectively for erg/s, MeV, erg/s, and .
The grey colored area denotes the parameter space where the transition to supercriticality occurs due to the nonlinear feedback loops leading to efficient pion production. As we show in §3, the selfconsistently produced gammaray spectrum starts resembling a typical GRB one, for or . Thus, our model is equivalent to the WB model in the sense that efficient pion production on a Bandlike gammaray spectrum is ensured for roughly similar parameter values, but with a fundamental difference: here the gammaray spectrum is not assumed a priori but is produced selfconsistently by a series of processes, which we describe in detail in §3.2.
3 Numerical investigation
3.1 Numerical code
In the previous section we showed by analytical means that the transition to supercriticality is ensured for a wide range of parameter values. To verify this we employ the timedependent numerical code as presented in Dimitrakoudis et al. (2012a) – hereafter DMPR12, that follows the evolution of protons, neutrons, secondary pairs, photons and neutrinos by solving the coupled differential equations that describe the various distributions. The coupling of energy losses and injection introduces a selfconsistency in this approach that allows the study of the system at various conditions, e.g. in the presence of nonlinear electromagnetic (EM) cascades and other feedback loops (see also Petropoulou & Mastichiadis 2012b for a relevant discussion). As a word of caution we stress that the aforementioned loops can be fully understood only if the coupled kinetic equation approach is used. The often used Monte Carlo techniques are intrinsically linear and fail to capture complex, nonlinear effects such as this. While Monte Carlo codes are an excellent tool for the description of the system in the subcritical regime, should the choice of parameters drive the system into the supercritical regime the results of such codes may be in error.
We assume that protons are being injected in the source at a constant rate given by
(17) 
where , , , and is the time measured in the comoving frame in units. Protons, as well as secondary particles, are allowed to leave the emission region in an average time . This may account in an approximate way for the expansion of the source, since the steadystate particle distributions derived by solving a kinetic equation containing a physical escape term or an adiabatic loss term are similar. All particles, primary and secondary, lose energy through various processes. Although details can be found in DMPR12, for the sake of completeness, we summarize here the physical processes that are included in the code:

protonphoton pair production (photopair)

protonphoton pion production (photopion)

neutronphoton pion production

proton synchrotron radiation

pion, kaon, muon and electron synchrotron radiation

synchrotron selfabsorption

electron inverse Compton scattering

photonphoton pair production

electronpositron pair annihilation

Compton scattering of photons by cooled pairs
Photohadronic interactions are modelled using the results of Monte Carlo simulations. In particular, for BetheHeitler pair production the Monte Carlo results by Protheroe & Johnson (1996) were used (see also Mastichiadis et al. 2005). Photopion interactions were incorporated in the timedependent code by using the results of the Monte Carlo event generator SOPHIA (Mücke et al., 2000). Synchrotron radiation of charged pions and muons was not included in the version of the code presented in DMPR12 and for the exact treatment we refer the reader to Dimitrakoudis et al. (2014). Pairs that cool down to Lorentz factors contribute to the Thomson depth and they are treated as a separate population. Following Lightman & Zdziarski (1987), we assume that this population thermalizes at a temperature , where . For the pair annihilation and photon downscattering processes we followed Coppi & Blandford (1990) and Lightman & Zdziarski (1987), respectively (for more details see Mastichiadis & Kirk 1995).
The photon escape timescale , in particular, is modelled as
(18) 
where is the photon energy in units and is a function that equals unity for , it decreases as for and it becomes zero for . Moreover, and is the KleinNishina cross section. On the one hand, the above expression takes into account in an approximate way the fact that photons may be ‘trapped’ in the source for longer than one crossing time, i.e. for . On the other hand, eq. (18) does not take into account the effect of the expansion of the source during the photon escape. If the source expands on a dynamical time, the photon escape time is found to be (Giannios, 2006). Note, however, that in the examples shown here , making thus the effect of expansion rather modest.
3.2 Results
We investigated in total 10 parameter sets that were divided in two groups, namely A and B. All simulations in groups A are obtained for erg/s, , and s, whereas for group B s. In each group we performed 5 simulations with Lorentz factors varying between and with a logarithmic step of 0.1. All other parameters that are used as an input for the numerical code, i.e. , , , and , are then derived using eqs. (1), (3), (6)(7), and (5) respectively. These are summarized in Table 1. In all cases we let the system reach a steadystate, where the photon and neutrino emission as well as the proton and neutron energy distributions were then calculated. Although, in most cases, a steadystate is achieved in , a timedependent treatment of the GRB emission that is intrinsically variable, is more adequate and it will be the subject of a future work.
#  (G)  (cm)  ^{a}  
Group A: s  
1  
2  
3  
4  
5  
Group B: s  
6  
7  
8  
9  
10 

It is calculated using the first branch of eq. (13) without the logarithmic dependance.
3.2.1 Emission spectra
Figure 3 shows the observed multiwavelength photon spectra (solid lines) obtained from a single GRB pulse for erg/s, , and s (Group A) and for three indicative values of , i.e. (black lines), (blue lines) and (black lines). For comparison reasons, the total neutrino^{4}^{4}4We refer to the sum of electron/muon neutrino and antineutrino fluxes as the total neutrino flux. (dotted lines) and proton escaping (dashed lines) energy spectra are overplotted. Since is kept fixed, the proton injection compactness increases for decreasing (see also eq. (5)). In particular, it increases gradually from (lower curve) to (upper curve) with increaments of 0.5 in logarithm. The ratio calculated using eq. (14) for each case is given in label of Fig. 3.
Spectra shown with grey lines are obtained for , which for the particular choice of parameters, leads to low and . This is a typical example of emission signatures obtained when the system operates in the subcritical regime. The photon emission is characterized by the following: the MW spectra are dominated by the proton synchrotron component, i.e. synchrotron radiation is the dominant energy loss mechanism for UHE protons, and the radiative efficiency, which is defined as , is low in agreement with typical proton synchrotron emission models (see e.g. Mücke & Protheroe (2001)). In the subcritical regime, the proton spectra at steady state are the same as at injection because of negligible energy losses. Note that the sharp cutoff of the proton spectra at reflects the injection spectrum . The high energy cutoff would have been smoother if we were to use a more physically motivated injection spectrum, e.g. . The produced neutrinos in this regime are the result of photopion interaction of protons with their own emitted synchrotron radiation (see also DMPR12). However, the photopion energy loss rate of protons is very low and the neutrino emission is supressed. Note that the peak neutrino luminosity is .
As the injection compactness increases, at some point it will exceed the critical value (see eq. (13)) and the system will undergo a phase transition, which can be easily identified by a radical change of the spectral shape and an abrupt increase of the emitted luminosity (see black and blue lines in Fig. 3). Note that all spectra obtained in the supercritical regime are also characterized by in agreement with the analysis of §2.2. The underlying reason for this abrupt transition is the instability of spontaneous ray quenching that redistributes the energy from the ray energy band to lower energy parts of the spectrum through an EM cascade (Stawarz & Kirk, 2007; Petropoulou & Mastichiadis, 2011). Because of the increased production of low energy photons in the source, it is the photopion/photopair energy loss channels that are favoured over the the proton synchrotron one while in the supercritical regime. Thus, the EM cascade initiated by the instability of spontaneous ray quenching is further supported by the injection of secondary pairs, which are highly relativistic and produced by the photopair and photopion interactions of UHE protons with the automatically generated soft photons. In other words, once protonproduced rays reach a certain compactness they are spontanteously absorbed, giving rise to electronpositron pairs and radiation, which causes more proton cooling via photopair and photopion processes; and eventually more rays, thus sustaining the loop, which is illustrated in Fig. 4.
All the processes discussed above lead naturally to a large number density of cooled pairs () that is related to the source’s Thomson depth () through . Since the number of secondaries produced in the EM cascade is related to the proton injection compactness, we expect that larger values of lead to larger optical depths. This trend is exemplified in Fig. 5 where we plot the Thomson optical depth versus at the injection (circles) and at the steady state (squares) for the same parameters as in Fig. 3. The hatched area corresponds to the supercritical regime, while indicative values of (eq. 14) are marked on the plot for comparison reasons. At injection we assumed that any electrons present in the source are these related to the injected protons for conservation of neutrality and, hence, . At the steady state we find for , and this is another manifestation of the increased secondary pair injection while the system lies in the supercritical regime.
If the value of is such as to result in (blue lines in Fig. 3), then the process that determines the spectral shape is photon Comptonization by cooled electrons. In particular, our numerical simulations indicate that the peak energy of the photon spectrum scales approximately as , which is in good agreement with the dependance found by the solving the Kompaneets equation (Kompaneets, 1956). Although similar work can be found in the literature (e.g. Arons (1971); Illarionov & Syunyaev (1972); Lightman et al. (1981)), we present in Appendix B an analytic derivation for the case of continuous powerlaw photon injection that better describes the physical problem under investigation, and further supports our numerical results. Note that if we were to neglect the effects of photon downscattering, then all spectra obtained in the supercritical regime would have had the universal shape of a spectrum that peaks at GeV in the observer’s frame (see e.g. blue line in Fig. 3), which reminds of the powerlaw underlying component seen in several bright GRBs (e.g. 080319B, 090902B, 090926A) detected with Fermi (Racusin et al., 2008; Abdo et al., 2009; Ackermann et al., 2011). Cascade emission produced through proton interactions with the MeV photons of the GRB were suggested as an alternative explanation for this emission (e.g. Asano et al. 2010). Interestingly, we derive similar photon spectral shapes which are also the result of a cascade but with the key difference that a photon field is not required ab initio; the targets are provided by the quenching loop. In the succession of spectra shown in Fig. 3 only those obtained for correspond to large enough and to start resembling to a Bandlike photon spectrum (Band et al., 1993). The exact spectral shapes, however, should be considered with caution, since a better description for the cooling of pairs below is required and which we plan to address in the future.
As the system is driven deeper to the supercritical regime the photopion energy losses become gradually more significant and the energy drained from UHE protons is transfered to photons and secondary particles, such as pairs and neutrinos. The abrupt increase of the neutrino fluence is an additional sign of the transition to the supercritical regime – see grey and black lines in Fig. 3. As the proton injection compactness progressively increases (bottom to top) the neutrino spectrum becomes flatter in units, and it extends to lower energies tracing the evolution of the peak photon luminosity. The neutrino spectra for share many common features with those obtained in studies where the photon target field is modelled by an ab initio Band spectrum (e.g. Murase 2008; Baerwald et al. 2011, 2013; Petropoulou 2014).
Summarizing, conditions leading to high photon and neutrino luminosities cause unavoidable cooling of UHE protons (see also Asano 2005) that obtain a steeper powerlaw spectrum than that of injection, i.e. , where . Figure 6 compares the observed differential luminosity of escaping protons and neutrons for the same parameters as in Fig. 3. In all three cases, we find that protons with energies eV are unaffected by cooling due to photohadronic processes and actually serve as a large energy reservoir. However, their contribution to the spectral and temporal properties of the gammaray spectra is negligible both in the subcritical and supercritical regimes, and this can be understood as follows: protons with eV emit mainly through synchrotron at low photon energies, e.g. and with a much lower luminosity than the gammaray one^{5}^{5}5The gammaray luminosity in the subcritical regime is given by the peak luminosity of the proton synchrotron component, while in the supercritical regime is given by the peak of the cascade emission component (see Fig. 3).. In this context it is the high energy part of the proton distribution that is ‘active’. The luminosity carried by neutrons, which are produced via the photopion channel , increases as the cooling of UHE protons becomes progressively more significant, i.e. as the conditions lead the system deeper into the supercritical regime. At the same time, the peak of the neutron energy spectrum moves towards lower energies, similarly to the cooling break energy of the proton energy spectrum, and the neutron distribution can be described by the same powerlaw as cooled protons, i.e., with index . In the optically thin limit for photopion interactions the produced neutron distribution would follow the proton injection spectrum. Here, the steepening of the neutron spectrum is the result of efficient neutron cooling through before escape from the emission region. Contrary to protons, neutrons are not confined by magnetic fields and their decay time is large enough to allow them to escape freely before converting through decay into protons. These will propagate as UHECRs into the intergalactic medium having an injection spectrum similar to that of their parent population (Kirk & Mastichiadis, 1989; Begelman et al., 1990; Giovanoni & Kazanas, 1990; Mannheim et al., 2001; Atoyan & Dermer, 2003). Other possible escape mechanisms of UHECRs are discussed in Baerwald et al. (2013); Asano & Meszaros (2014).
Having explained the basic features of the photon, neutrino and proton spectra as the system is driven progressively from the subcritical to the supercritical regime for a particular parameter set (Group A), we proceed to investigate the role of other parameters, such as – see Fig. 7. If all other parameters are kept fixed, faster variability is translated to higher proton injection compactness (), higher magnetic field strength (), smaller emission region () but approximately constant , as the only dependance on comes through the logarithmic term (see eq. (13)). Thus, for the same but smaller the system is driven more deep to the supercritical regime (panel (a) in Fig. 7). The effect of the stronger magnetic field for s is also imprinted on the cutoff energy of the neutrino spectrum, which moves from PeV to 10 PeV (panel (b)). The gammaray emission produced for in both cases is shown in Fig. 8. The bowties show the range of observed values for the GRB photon indices, namely and (Preece et al., 2000), and are plotted for guiding the eye. Note the change of the gammaray spectrum below the peak, which becomes harder as .
3.2.2 Efficiency
A robust sign of the transition to supercriticality is the abrupt increase of the photon, neutron and neutrino luminosities. The efficient conversion of energy originally stored in relativistic protons into radiation can be the result of other underlying feedback loops, such as the ‘PPSloop’ (Kirk & Mastichiadis, 1992), which was also studied in the framework of GRB prompt emission (Mastichiadis & Kazanas, 2006, 2009). Here we define the efficiency of the th component as , where , (photons), (neutrinos), (neutrons), (protons) and (protons with energies eV). The efficiency as a function of the Lorentz factor is shown in Fig. 9 with panels (a) and (b) corresponding to Groups A and B, respectively. Few things are worth commenting:

The luminosity of the proton component is the dominant one. The photon luminosity becomes comparable to the luminosity carried by the proton component only for parameters that drive the system deep into the supercritical regime (see e.g. Fig. 1), where we typically find that or .

The abrupt increase of , and is seen for (see also eq. (15)) and it marks the transition to supercriticality.

In panel (b) the abrupt increase of the neutrino, neutron and photon luminosities is not evident, as our simulations do not extend above , which in this case is .

The steep decrease of ( orders of magnitude for a 0.3 order of magnitude change in ) indicates the significant cooling of UHE protons.

For we find . For on the other hand, decreases faster than , since photopion interactions are the sole source of neutrinos contrary to photons, which are produced mainly via synchrotron radiation in the subcritical regime.

The efficiencies in neutrons and neutrinos are, generally, of the same order of magnitude. They show, however, different trends: remains approximately constant (panel (a)) or decreases (panel (b)) for smaller values of , whereas increases. This implies that neutrons interact with photons before escaping from the source and contribute to the neutrino production through the process .
3.2.3 Maximum proton energy
We continue our study on the emission features while the system is in the supercritical regime by considering a fiducial case where protons are not accelerated up to UHE, i.e. the maximum proton energy in the comoving frame is eV. One could imagine a scenario where the magnetization of the burst is very low (e.g. ) or the size of the emission region is small enough to confine higher energy protons.
To exemplify the above we adopt the following parameters: , erg/s, , s, and two values of the proton injection luminosity, erg/s and erg/s, which correspond to the subcritical and supercritical regime, respectively. Other parameters used are: G and cm.
The results for the subcritical and supercritical cases are summarized in panels (a) and (b) of Fig. 10, respectively. In the subcritical regime the various components of the overall photon spectrum can be identified, in contrast to the supercritical regime where the formation of the EM cascade blurs the emission signatures from individual processes. In panel (a) we plot the various contributions in order to demonstrate that in the TeV energy band (grey colored region) the gammaray emission is no longer dominated by the proton synchrotron component (for comparison see grey line in Fig. 3). Instead, it is the synchrotron radiation from secondary pairs produced through the BetheHeitler and photopion processes that is being emitted as Tev gammarays and that it is going to initiate the instability of automatic photon quenching (for relevant discussion see also Petropoulou & Mastichiadis 2012b).
Note that the emission signatures of hadronic plasmas in the subcritical regime may differ significantly for various parameter sets, as we exemplified in Figs. 3 and 10. However, the photon and neutrino emission produced when the system is driven to the supercritical regime, and in particular for parameters which ensure , is rather ‘universal’. Besides the overall energetics, the photon SEDs shown in Figs. 3 (black line), 8 and 10 have similar features.
We also found that for lower values of , a higher value is generally required to enter the supercritical regime. Given that in this regime approximetaly goes to gammarays and neutrinos (see also Fig. 9), the above conditions lead inevitably to bright photon and neutrino bursts. For the case shown in Fig. 10 for example, the gammaray fluence emitted in the supercritical regime would be erg/cm, if we were to assume a duration of s. This would place such an event at the highfluence tail of the Fermi distribution (von Kienlin et al., 2014).
4 Discussion
Hadronic models have served for a long time as alternatives to the more popular leptonic ones for AGN and GRB highenergy emission. One of their unique features is the prediction of copious neutrino emission which can be produced alongside the photon spectrum and is particularly attractive because it allows for hadronic models to be further tested. These suffer, however, from low efficiencies, as hadrons have in general long cooling timescales.
In earlier works (Mastichiadis et al., 2005; Dimitrakoudis et al., 2012b; Petropoulou & Mastichiadis, 2012a, b) that were targeting AGN highenergy emission, we showed that we can divide the parameter space of hadronic plasmas into two regimes. In the first one, which we shall call ‘subcritical’, protons carry the majority of the energy while a small amount is radiated away mainly by proton synchrotron radiation; photon induced processes like photopair and photopion carry even less luminosity – this is especially true if one assumes that there are no ambient photons illuminating the source. In this regime the system is inefficient, i.e. protons lose a very small part of their energy to secondaries. The other regime, the ‘supercritical’ one, is separated by a sharp boundary in phase space from the subcritical one and is characterised by exactly the opposite trend. Here, processes like photopair and photopion production dominate the losses and essentially drain the protons of their stored energy giving it to secondaries, thus increasing the efficiency to high values. Because of the radical change not only in the efficiency but also in the emission signatures of such plasmas, the transition from the subto the supercritical regime can be characterised as a ‘phase’ transition, with the underlying reason being the existence of nonlinear feeback loops.
In the present study we reexamined the radiative signatures of hadronic plasmas in the context of GRBs. The key question we wanted to answer was whether the energetics required for the supercriticality are compatible with typical GRB parameters. For this, we adopted the usual GRB hadronic picture used in the literature, i.e. we assumed the injection of high energy protons having a power law distribution in a source of a given size that contains a certain magnetic field and adopted parameter values relevant to GRB sources. However, departing from the usual GRB assumptions, neither external photons nor an extra population of accelerated relativistic electrons were considered. This choice minimized the number of free parameters to five: , , , and ; all other quantities, such as the proton injection luminosity () and the maximum energy () of their distribution, can be expressed in terms of these parameters (see eqs. (1), (3), (5), (6), (7)). For the derivation of we assumed that proton acceleration is fast, i.e. , where is the proton gyration timescale. For example, stochastic Fermi acceleration due to magnetic turbulence (e.g. Waxman 1995a; Dermer & Humi 2001) as well as relativistic magnetic reconnection (e.g. Giannios 2010) have been suggested as viable processes for fast, UHE proton acceleration. Whether or not the accelerated protons can form a powerlaw distribution remains to be shown, although there are some indications for the formation of a highenergy proton powerlaw tail in simulations of relativistic reconnection in electronion plasmas (private communication with Dr. L. Sironi). Note that the same kind of numerical studies performed for pair plasmas clearly show the formation of a powerlaw electron distribution (Sironi & Spitkovsky, 2014).
As a tool for our study we employed a recently developed numerical code (Dimitrakoudis et al., 2012a) for solving the system of spatially averaged kinetic equations, that describes the coupling between protons and their stable byproducts, namely photons, electrons (and positrons), neutrons and neutrinos. This, in contrast to most of the work performed on hadronic plasmas thus far, has allowed us to make a selfconsistent study of the evolution of the system by keeping track of the energy lost and gained by the various species. We have confirmed that this coupling between the species is the key to understanding its behaviour as only with this scheme can one follow the development and growth of nonlinear loops that eventually lead the system to supercriticality.
Altough there are various loops that may cause the phase transition from sub to supercriticality (e.g. Kirk & Mastichiadis 1992; Kazanas et al. 2002; Mastichiadis et al. 2005), our analysis presented in §2 and in Appendix A shows that ray quenching (Stawarz & Kirk, 2007; Petropoulou & Mastichiadis, 2011) is the leading one for the parameters used. This means that once protonproduced rays reach a certain compactness they are spontanteously absorbed, giving rise to electronpositron pairs and radiation, which causes more proton cooling via photopair and photopion processes; and eventually more rays, thus sustaining the loop. The cycle will continue until protons are drained of their energy and this, depending on the initial conditions, can lead the system either to steady state or to a limit cycle behaviour – in the latter case energy is gradually built into protons and is abruptly released in a few crossing times. Note that all the examples presented in §3 were obtained for parameters that led quickly (in 12 dynamical times) to a steady state.
As we have shown in §2, it is the existence of a critical ray compactness that makes the proton compactness the most important parameter in our study. On the one hand, if we use a similar definition to the radiation compactness and relate the observed proton luminosity, which is assumed to be a fraction of the jet’s kinetic luminosity , to the one measured in the comoving frame of the flow through a Lorentz transformation, we find . On the other hand, we showed that the critical ray compactness is translated to a critical value for the proton one, which roughly scales as . The combination of the and dependances favours supercriticality for most parameter values, except for flows with high bulk Lorentz factors. These results are summarized in Fig. 1, which answers in the most satisfactory manner our earlier posed key question.
Apart from the above, there are some far reaching consequences of this model when applied to GRBs. First, the efficiency of photons and neutrinos becomes quite high, reaching 0.10.4 of the total available proton luminosity. Moreover, the high number of electronpositron pairs created as secondaries cool to low energies producing a high Thomson optical depth (ranging from a few up to ) which can downscatter high energy photons producing a bump at in the rest frame of the flow (see Appendix B). It is noteworthy that a similar idea for explaining the peak of the GRB emission was proposed by Brainerd (1994), although in a different context. Moreover, the boundary in the phase space between the two regimes is very sharp, in the sense that a small perturbation in one of the proton injection parameters, while the system is still in the subcritical regime, can push it over to the supercritical one. The transition is in most cases very abrupt and it manifests itself with a photon flare which lasts a few crossing times (Dimitrakoudis et al., 2012a; Petropoulou & Mastichiadis, 2012b).
Clearly there are open questions that one needs to address before the present model can successfully explain the GRB phenomenology, such as gammaray spectra and timevariability, which we plan to investigate in a forthcoming publication. For the former, we plan to include an additional equation for lowenergy electrons () and follow the cooling/heating due to Compton process in more detail. Regarding the latter, one of the central issues that has to be addressed is whether or not our hadronic model can produce the fast variability observed in rays. Preliminary results of photon lightcurves derived in the supercritical regime assuming a variable proton injection rate show that rapid variability (of the order of the source crossing time) can be obtained. Note that even for a constant proton injection rate there are regimes in the parameter space that are relevant to GRBs and lead to a limit cycle behaviour (Petropoulou & Mastichiadis, 2012b). Thus, in cases where the proton injection is variable we find that the structure of the gammaray light curves is complex due to the superposition of the intrinsic periodicity of the hadronic plasma and of the variability pattern of the ‘external’ source of proton injection.
5 Summary
We showed that the injection of highenergy protons, e.g. eV, with luminosities erg/s in a region which is part of a GRBlike flow with bulk Lorentz factor can lead to an abrupt energy transfer from protons to photons and neutrinos, which may carry approximately of the injected proton luminosity. We also found that the proton injection luminosity that is required for triggering the efficient cooling of protons is higher, e.g. reaching erg/s, for proton distributions extending to lower energies, e.g. . In this case, our model has the testable prediction of a contemporaneous bright burst in rays and highenergy neutrinos. If the neutrino nondetection from Fermi bright GRBs with IceCube will be established (e.g. Abbasi et al. 2012; He et al. 2012; Liu & Wang 2013), in the context of our model, it will mean that the conditions in sources where proton acceleration to UHE is not possible, should be such as not to drive the system to supercriticality, e.g. erg/s.
We showed that in our framework the gammaray photon spectrum is selfconsistently determined. Its shape is sensitive on how deep in the supercritical regime the system is driven. When the proton luminosity is marginally above the critical value, the photon spectrum peaks at GeV, whereas for higher proton luminosities it is modified due to Comptonization and appears more as a Band spectrum. Not only different values of the proton luminosity but also small changes in one of the other parameters, and in particular of the bulk Lorentz factor, may lead to photon spectra ranging between the two shapes described above.
Summarizing, we showed that supercriticalities are a generic feature of hadronic models, as they manifest themselves for a wide range of parameters; from those relevant to AGN highenergy emission to those relevant to GRBs. In the context of the latter, such supercriticalities can lead naturally to bursts of gammarays that share several properties with the typical GRB prompt emission, such as the observed photon luminosity, the Bandlike spectra and the time duration. These offer also a unique way of transferring energy from protons to photons and neutrinos in a very efficient way, and of overcoming the usual lowefficiency problem of hadronic models.
Acknowledgements
Support for this work was provided by NASA through Einstein Postdoctoral Fellowship grant number PF 140113 awarded by the Chandra Xray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS803060. DG acknowledges support from the Fermi 6 cycle grant number 61122.
References
 Aartsen et al. (2014a) Aartsen M. G. et al., 2014a, PhysRevD, 89, 102001
 Aartsen et al. (2014b) Aartsen M. G. et al., 2014b, ArXiv eprints
 Abbasi et al. (2012) Abbasi R. et al., 2012, Nature, 484, 351
 Abbasi et al. (2010) Abbasi R. et al., 2010, ApJ, 710, 346
 Abdo et al. (2009) Abdo A. A. et al., 2009, ApJL, 706, L138
 Ackermann et al. (2011) Ackermann M. et al., 2011, ApJ, 729, 114
 Arons (1971) Arons J., 1971, ApJ, 164, 437
 Asano (2005) Asano K., 2005, ApJ, 623, 967
 Asano et al. (2009) Asano K., Guiriec S., Mészáros P., 2009, ApJL, 705, L191
 Asano & Inoue (2007) Asano K., Inoue S., 2007, ApJ, 671, 645
 Asano et al. (2010) Asano K., Inoue S., Mészáros P., 2010, ApJL, 725, L121
 Asano & Mészáros (2013) Asano K., Mészáros P., 2013, JCAP, 9, 8
 Asano & Meszaros (2014) Asano K., Meszaros P., 2014, ArXiv eprints
 Atoyan & Dermer (2003) Atoyan A. M., Dermer C. D., 2003, ApJ, 586, 79
 Baerwald et al. (2013) Baerwald P., Bustamante M., Winter W., 2013, ApJ, 768, 186
 Baerwald et al. (2014) Baerwald P., Bustamante M., Winter W., 2014, ArXiv eprints
 Baerwald et al. (2011) Baerwald P., Hümmer S., Winter W., 2011, PhysRevD, 83, 067303
 Band et al. (1993) Band D. et al., 1993, ApJ, 413, 281
 Band et al. (2009) Band D. L. et al., 2009, ApJ, 701, 1673
 Begelman et al. (1990) Begelman M. C., Rudak B., Sikora M., 1990, ApJ, 362, 38
 Beloborodov (2010) Beloborodov A. M., 2010, MNRAS, 407, 1033
 Bloom et al. (2003) Bloom J. S., Frail D. A., Kulkarni S. R., 2003, ApJ, 594, 674
 Boettcher (2010) Boettcher M., 2010, ArXiv eprints
 Böttcher (2007) Böttcher M., 2007, Ap&SS, 309, 95
 Brainerd (1994) Brainerd J. J., 1994, ApJ, 428, 21
 Chiang & Dermer (1999) Chiang J., Dermer C. D., 1999, ApJ, 512, 699
 Coppi & Blandford (1990) Coppi P. S., Blandford R. D., 1990, MNRAS, 245, 453
 Crider et al. (1997) Crider A. et al., 1997, ApJL, 479, L39
 Daigne et al. (2011) Daigne F., Bošnjak Ž., Dubus G., 2011, A&A, 526, A110
 Derishev et al. (2001) Derishev E. V., Kocharovsky V. V., Kocharovsky V. V., 2001, A&A, 372, 1071
 Dermer & Atoyan (2003) Dermer C. D., Atoyan A., 2003, Physical Review Letters, 91, 071102
 Dermer & Atoyan (2006) Dermer C. D., Atoyan A., 2006, New Journal of Physics, 8, 122
 Dermer & Humi (2001) Dermer C. D., Humi M., 2001, ApJ, 556, 479
 Dimitrakoudis et al. (2012a) Dimitrakoudis S., Mastichiadis A., Protheroe R. J., Reimer A., 2012a, A&A, 546, A120
 Dimitrakoudis et al. (2012b) Dimitrakoudis S., Petropoulou M., Mastichiadis A., 2012b, International Journal of Modern Physics Conference Series, 8, 19
 Dimitrakoudis et al. (2014) Dimitrakoudis S., Petropoulou M., Mastichiadis A., 2014, Astroparticle Physics, 54, 61
 Gao et al. (2012) Gao S., Asano K., Mészáros P., 2012, JCAP, 11, 58
 Giannios (2006) Giannios D., 2006, A&A, 457, 763
 Giannios (2010) Giannios D., 2010, MNRAS, 408, L46
 Giannios (2012) Giannios D., 2012, MNRAS, 422, 3092
 Giovanoni & Kazanas (1990) Giovanoni P. M., Kazanas D., 1990, Nature, 345, 319
 Goldstein et al. (2012) Goldstein A. et al., 2012, ApJS, 199, 19
 Goodman (1986) Goodman J., 1986, ApJL, 308, L47
 Guetta et al. (2004) Guetta D., Hooper D., AlvarezMun~Iz J., Halzen F., Reuveni E., 2004, Astroparticle Physics, 20, 429
 He et al. (2012) He H.N., Liu R.Y., Wang X.Y., Nagataki S., Murase K., Dai Z.G., 2012, ApJ, 752, 29
 Hillas (1984) Hillas A. M., 1984, ARA&A, 22, 425
 IceCube Collaboration (2013) IceCube Collaboration, 2013, Science, 342
 Illarionov & Syunyaev (1972) Illarionov A. F., Syunyaev R. A., 1972, SvA, 16, 45
 Kardashev (1962) Kardashev N. S., 1962, SvA, 6, 317
 Katz (1994) Katz J. I., 1994, ApJL, 432, L107
 Kazanas et al. (2002) Kazanas D., Georganopoulos M., Mastichiadis A., 2002, ApJL, 578, L15
 Kirk & Mastichiadis (1989) Kirk J. G., Mastichiadis A., 1989, A&A, 213, 75
 Kirk & Mastichiadis (1992) Kirk J. G., Mastichiadis A., 1992, Nature, 360, 135
 Kocevski & Butler (2008) Kocevski D., Butler N., 2008, ApJ, 680, 531
 Kompaneets (1956) Kompaneets A. S., 1956, Zh.E.F.T, 31, 876
 Lazzati et al. (2000) Lazzati D., Ghisellini G., Celotti A., Rees M. J., 2000, ApJL, 529, L17
 Lightman et al. (1981) Lightman A. P., Lamb D. Q., Rybicki G. B., 1981, ApJ, 248, 738
 Lightman & Zdziarski (1987) Lightman A. P., Zdziarski A. A., 1987, ApJ, 319, 643
 Liu & Wang (2013) Liu R.Y., Wang X.Y., 2013, ApJ, 766, 73
 Lloyd & Petrosian (2000) Lloyd N. M., Petrosian V., 2000, ApJ, 543, 722
 Mannheim et al. (2001) Mannheim K., Protheroe R. J., Rachen J. P., 2001, PhysRevD, 63, 023003
 Mastichiadis & Kazanas (2006) Mastichiadis A., Kazanas D., 2006, ApJ, 645, 416
 Mastichiadis & Kazanas (2009) Mastichiadis A., Kazanas D., 2009, ApJL, 694, L54
 Mastichiadis & Kirk (1995) Mastichiadis A., Kirk J. G., 1995, A&A, 295, 613
 Mastichiadis et al. (2005) Mastichiadis A., Protheroe R. J., Kirk J. G., 2005, A&A, 433, 765
 Medvedev (2000) Medvedev M. V., 2000, ApJ, 540, 704
 Mészáros & Rees (2000) Mészáros P., Rees M. J., 2000, ApJ, 530, 292
 Mücke et al. (2000) Mücke A., Engel R., Rachen J. P., Protheroe R. J., Stanev T., 2000, Computer Physics Communications, 124, 290
 Mücke & Protheroe (2001) Mücke A., Protheroe R. J., 2001, Astroparticle Physics, 15, 121
 Murase (2008) Murase K., 2008, PhysRevD, 78, 101302
 Murase et al. (2012) Murase K., Asano K., Terasawa T., Mészáros P., 2012, ApJ, 746, 164
 Murase et al. (2008) Murase K., Ioka K., Nagataki S., Nakamura T., 2008, PhysRevD, 78, 023005
 Nakar & Piran (2002) Nakar E., Piran T., 2002, MNRAS, 331, 40
 Norris et al. (1996) Norris J. P., Nemiroff R. J., Bonnell J. T., Scargle J. D., Kouveliotou C., Paciesas W. S., Meegan C. A., Fishman G. J., 1996, ApJ, 459, 393
 Panaitescu & Mészáros (2000) Panaitescu A., Mészáros P., 2000, ApJL, 544, L17
 Pe’er et al. (2006) Pe’er A., Mészáros P., Rees M. J., 2006, ApJ, 642, 995
 Pe’er & Zhang (2006) Pe’er A., Zhang B., 2006, ApJ, 653, 454
 Petropoulou (2014) Petropoulou M., 2014, submitted to MNRAS
 Petropoulou et al. (2014) Petropoulou M., Giannios D., Dimitrakoudis S., 2014, ArXiv eprints
 Petropoulou & Mastichiadis (2011) Petropoulou M., Mastichiadis A., 2011, A&A, 532, A11
 Petropoulou & Mastichiadis (2012a) Petropoulou M., Mastichiadis A., 2012a, MNRAS, 426, 462
 Petropoulou & Mastichiadis (2012b) Petropoulou M., Mastichiadis A., 2012b, MNRAS, 421, 2325
 Piran (2004) Piran T., 2004, Reviews of Modern Physics, 76, 1143
 Preece et al. (1998) Preece R. D., Briggs M. S., Mallozzi R. S., Pendleton G. N., Paciesas W. S., Band D. L., 1998, ApJL, 506, L23
 Preece et al. (2000) Preece R. D., Briggs M. S., Mallozzi R. S., Pendleton G. N., Paciesas W. S., Band D. L., 2000, ApJS, 126, 19
 Protheroe & Johnson (1996) Protheroe R. J., Johnson P. A., 1996, Astroparticle Physics, 4, 253
 Racusin et al. (2008) Racusin J. L. et al., 2008, Nature, 455, 183
 Reynoso (2014) Reynoso M. M., 2014, ArXiv eprints
 Sari et al. (1996) Sari R., Narayan R., Piran T., 1996, ApJ, 473, 204
 Sironi & Spitkovsky (2014) Sironi L., Spitkovsky A., 2014, ApJL, 783, L21
 Stawarz & Kirk (2007) Stawarz Ł., Kirk J. G., 2007, ApJL, 661, L17
 Stern & Svensson (1991) Stern B., Svensson R., 1991, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 391, Relativistic Hadrons in Cosmic Compact Objects, Zdziarski A. A., Sikora M., eds., p. 41
 Stern & Poutanen (2004) Stern B. E., Poutanen J., 2004, MNRAS, 352, L35
 Sultana et al. (2013) Sultana J., Kazanas D., Mastichiadis A., 2013, ApJ, 779, 16
 Tavani (1996) Tavani M., 1996, ApJ, 466, 768
 Thompson (1994) Thompson C., 1994, MNRAS, 270, 480
 Totani (1998) Totani T., 1998, ApJL, 509, L81
 Vietri (1995) Vietri M., 1995, ApJ, 453, 883
 Vietri (1997) Vietri M., 1997, Physical Review Letters, 78, 4328
 von Kienlin et al. (2014) von Kienlin A. et al., 2014, ApJS, 211, 13
 Wang et al. (2009) Wang X.Y., Li Z., Dai Z.G., Mészáros P., 2009, ApJL, 698, L98
 Waxman (1995a) Waxman E., 1995a, Physical Review Letters, 75, 386
 Waxman (1995b) Waxman E., 1995b, ApJL, 452, L1
 Waxman & Bahcall (1997) Waxman E., Bahcall J., 1997, Physical Review Letters, 78, 2292
 Zhang & Kumar (2013) Zhang B., Kumar P., 2013, Physical Review Letters, 110, 121101
 Zhang & Mészáros (2004) Zhang B., Mészáros P., 2004, International Journal of Modern Physics A, 19, 2385
Appendix A Derivation of critical proton compactness
Here we derive the expression for the critical proton compactness () in the case where the gammarays that are spontaneously absorbed are the result of proton synchrotron radiation. We derive the expression for a powerlaw proton distribution with slope . A similar analysis can be followed for different powerlaw indices.
The peak of the proton synchrotron spectrum in units appears at
(19) 
where , G, and . Here is the Lorentz factor of protons that cool within the dynamical timescale and it is given by eq. (9). The maximum Lorentz factor of protons is , where and are given by eqs. (6) and (7), respectively.

The condition is satisfied as long as
(22) Using eqs. (9) and (6) we find that only if
(23) Inspection of conditions (22) and (23) reveals that in the regime where the highenergy part of the proton distribution is affected by synchrotron losses, at least for most parameter values. For this, we set as long as condition (22) is satisfied.
The peak luminosity of the proton synchrotron spectrum in the comoving frame is then given by
(24) 
where
(25) 
is the total number of protons which is related to the proton injection rate as
(26) 
where , and . The above expression is exact only if . However, we find that is still a good approximation for and .
Using eq. (26) the proton injection compactness defined by eq. (5) can be written in terms of as
(27) 
Assuming that , where is the integrated ray luminosity, and using eqs. (24) and (27), the ray compactness is written as
(28) 
As shown by Petropoulou & Mastichiadis (2011), the critical ray compactness () is a function of the ray photon’s energy (see eq. (34) therein) having a minimum at the energy (in units)
(29) 
and increasing as for . The minimum value of is found to be . For , the critical ray compactness may be written as
(30) 
In general, the peak of the proton synchrotron spectrum () appears at . The transition to supercriticality occurs if . By combining eqs. (3), (19), (28), and (30) we find equivalently that
(31) 
Appendix B Peak energy  Thomson depth relation
The complete form of the Kompaneets equation (Kompaneets, 1956) is
(32) 
where , , and is the differential photon number density which is related to the photon occupation number as . In most astrophysically related cases and the induced emission term () in the Kompaneets equation may be safely neglected. Here we consider cases where and , i.e. the Kompaneets equation describes the ‘recoil effect’, where a photon loses energy through multiple scatterings with cold electrons.
By adding a source and an escape term, the Kompaneets equation now reads
(33) 
where is the time in units and is an approximate photon escape timescale given by
(34) 
and is given by eq. (21b) in Lightman & Zdziarski (1987). For our purposes, however, it is sufficient to use . Finally, the source term appearing in eq. (33) is given by
(35) 
We note that the Thomson depth is, in principal, a timedependent quantity. Using knowledge gained from the numerical study of the problem, according to which scales as with in the supercritical regime, we can consider to be constant.
Equation (33) is solved using the method of characteristics that transforms a partial differential equation (DE) into an ordinary DE along the characteristic curves or surfaces in two or three dimensional problems, respectively. In our case, the characteristic curve is given by
(36) 
same as in the case of synchrotron or/and inverse Compton (in the Thomoson regime) cooling of relativistic electrons (see e.g. Kardashev 1962). First, we find the solution to the homogeneous equation by setting . Along the characteristic curve of eq. (36), the PDE now reads