INVESTIGATING BROADBAND VARIABILITY OF THE TeV BLAZAR 1ES 1959+650

VERITAS is an array of four imaging atmospheric Cherenkov telescopes in southern Arizona, each with a 35 field of view. The array is sensitive to photons with energies from ∼100 GeV to more than 10 TeV and can detect a 1% Crab Nebula flux source at 5 standard deviations (σ) in less than 28 hr. The telescope array uses 12 m reflectors to focus dim, blue/UV Cherenkov light from gamma-ray and cosmic-ray interactions in the atmosphere onto cameras composed of 499 photomultiplier tubes (PMTs). More details on the VERITAS instrument can be found in Holder et al. (2006) and Weekes et al. (2002).

A historically high optical state of 1ES 1959+650 (R-band of 13.8 magnitude, as measured by the Super-LOTIS robotic telescope; Williams et al. 2008), prompted near-nightly VERITAS exposures over two dark runs³⁸ between 2012 April 17 and 2012 June 1 (MJD 56034-56079). Observations were carried out in wobble mode, with exposures taken at 05 offset in each of the four cardinal directions from 1ES 1959+650, in order to facilitate simultaneous background measurements (Fomin et al. 1994; Berge et al. 2007). The total exposure over the two dark runs resulted in 8.7 hr of quality-selected live time which was collected at an average zenith angle of 37°.

Elliptical moments of the recorded images are calculated and used to discriminate background cosmic-ray events from gamma-ray events. The data are first cleaned with "quality cuts," discarding any telescope images involving fewer than five PMTs or images with centroids at greater than 143 from the camera center. Additionally, each image is required to have a total "size" (a measure of total Cherenkov light collected by the camera) of more than ∼80 photoelectrons.

Single-telescope image widths and lengths are combined into mean-scaled-width (MSW) and mean-scaled-length (MSL) parameters for each array event, as described in Cogan (2008) and Daniel (2008). Only array events with 0.05 < MSW <1.25, 0.05 < MSL <1.1, a reconstructed height of shower maximum greater than 7 km above the array, and having a reconstruction direction within 01 of 1ES 1959+650 are kept as candidate signal (ON) events. The background (OFF) events are those that pass all aforementioned cuts and fall within 01 radius circular regions at the same radial distance as the source from the center of the camera. The source significance is calculated from the number of events falling in these ON and OFF regions according to Equation (17) of Li & Ma (1983). The analysis of the VHE source signal is confirmed with two independent analysis packages.

The total VERITAS exposure between 2012 April 17 and 2012 June 1 (MJD 56034 to 56079) resulted in 517 ON events and 1175 OFF events (corresponding to 410 excess events with an averaged background normalization α parameter of 0.0909) and an overall detection significance of 31σ. The VERITAS spectral data are derived with systematically coarser binning with increasing energy and are fit with a differential power law of the form  dN/ dE = N₀ × (E/E₀)^−Γ, where E₀ is fixed at 1 TeV and N₀ is the normalization parameter. Variability was detected during these observations, as can be seen by the source spectra in Figure 1 and in the top panel of the broadband light curve presented in Figure 2. The upper limits represent 95% confidence upper limits, calculated according to Rolke et al. (2005). A χ² test shows less than 6.4 × 10⁻¹² probability of a steady VHE flux. Table 1 contains a summary of the VHE analysis and spectral states of 1ES 1959+650 during the two dark runs (MJD 56034-56040 and 56064-56079) as well as, separately, the night where an elevated VHE state was detected (MJD 56067). This flare is excluded from the dark run analysis results. The results are shown with statistical errors only in Table 1. The systematic error on the energy scale is estimated between 20% and 35%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The hour-scale exposures during the first dark run show the blazar to be at an average flux of (8.6 ± 3.6) × 10⁻⁸ photons m⁻² s⁻¹ above the observational energy threshold of 315 GeV (approximately 8% of the Crab Nebula flux above this same threshold). The first dark-run observations are paired with two contemporaneous Swift exposures, described in Section 2.3.

On MJD 56067, VERITAS detected a short-lived VHE flare of 1ES 1959+650. These observations show the blazar flux to rise from ∼50% to 120% of the Crab Nebula in less than 30 minutes (see the top panel of Figure 3). The rise in flux was immediately followed by a decay, dropping back to ∼40% of the Crab Nebula approximately 90 minutes after the start of the event. The probability that the blazar flux on MJD 56067 was constant is less than 0.2% (see Figure 3; χ² = 29.5 with 11 degrees of freedom). The first 0.7 ks (12 minutes) of the VERITAS observations of the flaring event on MJD 56067 were matched with simultaneous Swift observations, described in Section 2.3. VERITAS continued to observe 1ES 1959+650 through 2012 June 1, detecting the source at an average of 12% Crab Nebula flux (1.5 ± 0.2 × 10⁻⁷ photons m⁻² s⁻¹ above 315 GeV).

Zoom In Zoom Out Reset image size

The Fermi Large Area Telescope (LAT) is a space-based telescope that typically monitors the entire high-energy gamma-ray sky from below 30 MeV to ∼300 GeV every three hours (Atwood et al. 2009). The instrument has better than 10% energy resolution, with an angular resolution of better than 015 for energies greater than 10 GeV.

Spectral analysis was completed for the period between (MJD 56054 to 56082) with the intention to search for an elevated gamma-ray state occurring contemporaneously with the elevated VHE states observed on MJD 56062 and MJD 56067. All analysis was completed with FermiTools v9r27p1. Events were extracted from a 30° radius region centered on the 1ES 1959+650 coordinates. "Diffuse class" events with zenith angles of <100° and energies between 300 MeV and 300 GeV were selected. In order to reduce contamination from Earth-limb gamma rays, only data taken while the rocking angle of the satellite was less than 52° were used. The significance and spectral parameters were calculated using the unbinned maximum-likelihood method gtlike with the P7SOURCE_V6 instrument-response functions. The background model was constructed to include nearby (<30° away) gamma-ray sources from the second Fermi LAT catalog (2FGL; Nolan et al. 2012) as well as diffuse emission.

As in the 2FGL catalog, a log-parabolic function was used for nearby sources with significant spectral curvature and a power law for those sources without spectral curvature. The spectral parameters of sources within 7° of 1ES 1959+650 were left free during fitting, while those outside of this range were held fixed to the 2FGL catalog values. The Galactic diffuse emission was modeled with the file gal_2yearp7v6_v0.fits and the isotropic emission component was modeled with the file iso_p7v6source.txt.³⁹

The analysis of LAT data between MJD 56054 and 56082 results in a test statistic (TS; Mattox et al. 1996) of 97 above 100 MeV. The spectral fitting shows the source to be in a slightly elevated state as compared to the 2FGL value (less than double the 2FGL integral flux of F_{1 − 100 GeV} = (8.8 ± 0.3) × 10⁻⁹ ph cm⁻² s⁻¹), with an integral flux of (1.6 ± 0.8) × 10⁻⁸ ph cm⁻² s⁻¹ above 1 GeV and a spectral index of 1.9 ± 0.1.

The data were binned in time to search for evidence of an elevated gamma-ray flux (see Figure 4). No evidence for variability on weekly timescales is found; the data are consistent with a steady flux at the 99.75% confidence level (χ² = 0.648 for 3 degrees of freedom). The LAT observations do not provide sufficient statistics for a more detailed investigation of variability, i.e., on daily timescales, which result in time-bin TS values of less than 9. Upper limits at 95% confidence are derived for daily LAT exposures of 1ES 1959+650 and are displayed in Figure 4 by downward pointing arrows.

Zoom In Zoom Out Reset image size

The Swift X-ray Telescope (XRT) is a space-based grazing incidence Wolter I telescope that focuses X-rays between 0.2 keV and 10 keV onto a 110 cm² CCD (Gehrels et al. 2004). The Swift telescope took 16 windowed timing (WT) exposures of 1ES 1959+650 between 2012 April 19 and 2012 June 1 (MJD 56036–56079), each between 0.5 ks and 1.5 ks long. The exposure on MJD 56067 is strictly simultaneous with VERITAS observations.

The data were analyzed using the HEASoft package Version 6.12. Rectangular source regions of length and width 45 pixels and 8 pixels, respectively, were used. Similarly sized regions of nearby source-free sky were used to estimate the background. The exposures were grouped to ensure a minimum of 20 counts per bin and were fitted with an absorbed power-law model of the form $F(E)_{{\rm PL}}= Ke^{-N_{\rm H\,\scriptsize{I}}\sigma (E)} (E/1\,$keV)^−α, with a free neutral hydrogen column density parameter N_H i, as in Furniss et al. (2013), to allow for additional absorption of soft X-rays by intervening gas in the vicinity of the blazar or along the line of sight. The model also contains a fitted normalization parameter K and the non-Thompson, energy-dependent photoelectric cross section σ(E), as taken from Morrison & McCammon (1983).

The flux and indices derived from the Swift XRT observations of 1ES 1959+650 are shown in the second and third panels from the top in Figure 2 and are summarized in Table 2. The fitted N_H i was consistently ∼two times higher than the Galactic value of 1 × 10²¹cm⁻², as measured by the LAB survey (Kalberla et al. 2005). The integral 2–10 keV flux recorded by the XRT ranged between 4.2 × 10⁻¹¹ erg cm⁻² s⁻¹ and 12.9 × 10⁻¹¹ erg cm⁻² s⁻¹ with an average flux of 7.5× 10⁻¹¹ erg cm⁻² s⁻¹, with photon indices ranging from α = 2.5 to 3.1.

The X-ray emission displayed by 1ES 1959+650 is relatively steady in the first five exposures. The exposure on 2012 May 20 (MJD 56067) is the only strictly simultaneous Swift observation with VERITAS, overlapping with the first 0.7 ks (∼12 minutes) of VERITAS observations during the VHE flare from 1ES 1959+650. The steady 0.3–10 keV count rate observed by XRT (as shown in the bottom panel of Figure 3) shows that at least the first 12 minutes of the VHE flaring episode is not matched with a simultaneous elevated X-ray state.

The Swift-XRT observations on MJD 56074 show the blazar to be in an elevated state with a flux level of (12.9 ± 0.1) × 10⁻¹¹ erg cm⁻² s⁻¹, nearly twice the average of 7.5× 10⁻¹¹ erg cm⁻² s⁻¹. This high state is observed to drop to approximately half the average (4.0 × 10⁻¹¹ erg cm⁻² s⁻¹) in less than two days. No contemporaneous high state is observed in the VERITAS data, but no firm conclusions can be drawn due to the non-simultaneous nature of the exposures.

Swift UVOT. The Swift-XRT observations were supplemented with simultaneous UVOT exposures taken in the UVW1, UVM2, and UVW2 bands (Poole et al. 2008). The Ultraviolet and Optical Telescope (UVOT) photometry is performed using the HEASoft program uvotsource. The source region consists of a single circle with 5'' radius, while the background region consists of several 15'' radius circles of the nearby sky lacking visible sources. The results are corrected for reddening using E(B − V) coefficients from Schlegel et al. (1998), with the Galactic extinction coefficients applied according to Fitzpatrick (1999). The largest source of error derived for the intrinsic flux points is due to the uncertainty in the reddening coefficients E(B − V). The UVOT W1, W2, and M2 flux values derived from the observations are shown in Figure 2 (fourth panel from the top), and summarized in detail in Table 3. These exposures show relatively steady flux values, with a small decrease of ∼30% up to MJD 56074, the day where an elevated X-ray state was observed with XRT.

Super-LOTIS. The Super-LOTIS robotic 0.6 m telescope located on Kitt Peak in Arizona took R-band exposures of 1ES 1959+650 between MJD 56034 and MJD 56080. During each night, three individual frames were acquired with the standard Johnson–Cousins R-band filter. Each image was reduced using an analysis pipeline which, after subtracting the bias and the dark current, combines the flat-fielded frames in a single image for each night. Aperture photometry with a circular aperture of 15'' was performed for both the blazar and each of the seven reference stars detailed in the Landessternwarte Heidelberg–Königstuhl catalogue⁴⁰ with a circular aperture of 15''. This aperture is large enough to encompass all the light enclosed in the irregular point-spread function. The local sky level is computed in a circular annulus of inner/outer radius of 18''/25''. The final flux values for 1ES 1959+650 are calculated by applying the photometric zero-point derived for each night, comparing the instrumental magnitude of the reference stars to the known magnitudes in the R band.

The R-band monitoring data from Super-LOTIS are shown in the bottom panel of Figure 2 and are summarized in Table 4. The observations show a relatively steady optical magnitude of between 13.8 and 14.0, with a conservative photometry error estimate of  ±0.1 optical magnitude.

iTelescope. V-band and R-band exposures were taken by the iTelescope between MJD 56028 and 56080. iTelescope is a robotic telescope system located in Nerpio, Spain.⁴¹ The telescopes used are twins (T07 and T18), and are each of 431 mm (17 inch) aperture at f/6.8. They employ an SBIG STL-1100M CCD camera. The V filter is a standard Johnson–Cousins set. The R filter is not standard and requires a color correction, where the addition of approximately 0.040 optical magnitude transfers the non-standard filter magnitudes to the standard Johnson–Cousins R. The data were reduced with MIRA Pro Version 7.0. The reduction is with standard aperture (radius of 5'') photometry, using the same standard stars as were used for the Super-LOTIS data reduction.

The V-band and R-band data are shown in the bottom panel of Figure 2 and summarized in Table 4. These observations show elevated luminosity in both the V and R bands (by approximately 0.2 optical magnitude) on two of the days where Super-LOTIS also provides R-band measurements approximately 7.5 hr after the iTelescope exposures were taken (shown in bold in Table 4). Comparison with the contemporaneous Super-LOTIS R-band measurements suggests that the blazar exhibits a small level of intranight variability. This fast variability occurs on the same night as the X-ray flux is observed to be high.

The broadband observations summarized above show a VHE flare on MJD 56067, where no elevated X-ray state is observed simultaneously for the first 12 minutes of the VHE flaring event. Additional variability is observed over the full window of observation, including an X-ray flux increase and intranight X-ray variability on MJD 56075. The X-ray flux was observed to drop over the next two days, with no corresponding (non-simultaneous) change in VHE flux observed.

We apply a time-independent SSC model to the relatively low and elevated flux states of 1ES 1959+650 observed on MJD 56064 and 56067. Both of these days have sufficient multiwavelength coverage to provide a full view of the broadband spectral energy distribution. Since the VHE flux of 1ES 1959+650 observed on MJD 56064 is consistent with the average flux from the two dark runs (excluding the VHE state on MJD 56067), the VHE spectrum is represented by the spectrum averaged over the two dark runs. Moreover, the average LAT spectrum is used to represent both the low and high states since no significant variability was detected within the LAT energy band.

The SSC model applied to the data is described in detail in Böttcher et al. (2013). This model is an equilibrium SSC model with emission originating from a spherical region of relativistic electrons with radius R. This region moves down the jet with a Lorentz factor Γ. The jet is oriented such that the angle with respect to the line of sight is θ_obs. In order to minimize the number of free parameters, the modeling is completed with θ_obs = 1/Γ, for which Γ = D, where D represents the Doppler boosting factor.

Electrons are injected into the spherical region with a power-law distribution of Q(γ) = Q₀γ^−q between the low- and high-energy cut-offs, γ_{min, max}. The electron distribution spectral indices used for 1ES 1959+650 are q = 2.7–2.8, which can be produced under acceleration in relativistic oblique shocks (Summerlin & Baring 2012). In order to reach an equilibrium state, the model evaluates the steady state produced when considering particle injection, radiative cooling, and particle escape. The particle escape is characterized with an efficiency factor η, such that the escape timescale t_esc = η R/c, with η = 1000 for this work, setting up an equilibrium scenario with a relatively long escape timescale for the relativistic particles. The variability timescale t_var is determined by the light crossing timescale of the emitting region (t_var = δR/c). According to this SSC model, the particle distribution streams along the jet with a kinetic power L_e. Synchrotron emission results from the presence of a tangled magnetic field B, with a Poynting flux luminosity of L_B. The parameters L_e and L_B allow the calculation of the equipartition parameter L_B/L_e.

A reasonable description to the low state (MJD 56064) is achieved with the parameters summarized in Table 5, and displayed by the solid line in Figure 5. Starting from this low state representation, two possible scenarios are explored to describe the elevated VHE state observed on MJD 56067. In the first realization (Scenario I, dotted line in Figure 5), the gamma-ray peak was shifted to a slightly higher energy by increasing the low-energy cut-off (γ_min). In order to keep the synchrotron peak at the same energy, the magnetic field is lowered. The lower magnetic field results in a lower synchrotron power, requiring an increase in overall jet power (L_e) to achieve the same synchrotron luminosity. This emission scenario, however, over-shoots the gamma-ray flux unless the radius of the emission region R is increased, lowering the compactness of the emission zone. The result of these parameter changes are shown by the dotted line in Figure 5, and summarized in the second column for Table 5. Under this representation, the variability timescale changes from 8.2 hr in the low state to 31 hr in Scenario I. This timescale is longer than the observed variability timescale on MJD 56067, where the VHE flux was observed to increase from 0.6× 10⁻⁶ photons cm⁻² s⁻¹ to 1.4× 10⁻⁶ photons cm⁻² s⁻¹ above 315 GeV in less than one hour.

Zoom In Zoom Out Reset image size

The eight parameters (including L_e) that describe the SSC emission scenario are known to be degenerate. As a result, the parameter changes described in Scenario I are not the only changes which will account for the difference between the SEDs observed on MJD 56064 and MJD 56067. Alternatively, in addition to the changes to the γ_min and magnetic field, an increase of the Doppler factor instead of a change to the emission region size can provide a similar result (Scenario II, dashed line in Figure 5, third column of Table 5). With this scenario, the variability timescale is still relatively short (5.1 hr), in agreement with the fast flux variability observed in the VHE band on MJD 56067.

As seen in Figure 5, both of the elevated state SSC scenarios predict an approximate doubling of the high-energy gamma-ray flux. This type of variability is impossible to rule out without detection of 1ES 1959+650 by the Fermi LAT on day-long timescales. The daily 95% confidence level upper limit on MJD 56067 (Figure 4) is more than double the LAT flux for the entire period (also used to derive the spectrum shown in the SED), indicating that the SSC-inferred LAT flux in an elevated state is still consistent with the observations. Therefore, the broadband SEDs of 1ES 1959+650 on MJD 56064 and MJD 56067 can be represented by a SSC emission model, necessitating multiple parameter changes from the high states relative to the low states in order to produce an elevated VHE state with no change to the synchrotron peak.

In this section, we present a possible scenario to describe the VHE variability detected during the contemporaneous multiwavelength observations of 1ES 1959+650. The main emphasis of our discussion is to show that the scenario can explain the data. Our choice in model is motivated, in part, by the evidence for intervening gas within the blazar. We apply a similar reflected emission model to that which was used to describe the "orphan" flaring activity of 1ES 1959+650 (Böttcher 2005). This model follows X-ray emission from a newly ejected component (blob) in the jet as it is reflected off dilute gas and/or dust in the vicinity of the jet. The reflected emission then re-enters the jet before the blob, which is moving down the jet, and reaches the location of the reflector. The application of this model is notably distinct from that applied in Böttcher (2005), where, for this application, the incident flux on the cloud is integrated over the time it takes the blob to pass the reflecting cloud instead of taken from a single short-lived X-ray flaring period. This integration over the blob's travel is necessary due to the assumed parsec-scale proximity of the blob to the reflecting cloud.

We assume that at a distance R_m = 1 R_{m, pc} pc from the central engine, moderately dense clouds of gas and dust (hereafter referred to as the "mirror") intercept the synchrotron emission from portions of the jet located inside R_m, and reprocess part of this flux back into the jet trajectory. Following Böttcher (2005), the distance R_m can be related to the observed time delay between the emergence of a new jet component from the core, and the (observer's frame) time at which the new component intercepts the mirror. For a characteristic bulk Lorentz factor of the new component of Γ = 10 Γ₁, and a time delay of, for example, Δt ∼ 5 × 10⁵ s:

$$\begin{equation} \Delta t \sim {R_m \over 2 \, \Gamma ^2 \, c} \sim 5 \times 10^5 \, R_{m, pc} \, \Gamma _1^{-2} \; {\rm s}, \end{equation} \tag{ 1 }$$

according to which the new component would have emerged around MJD 56062 for Γ ∼ 10 and R_m ∼ 1 pc.

The accumulated and reprocessed jet–synchrotron flux will be intercepted by the blob within the time interval between emitting the first photons at the time of emergence of the new component, and intercepting the location of the cloud. This time is also characteristic of the time during which the cloud receives this flux, and can be estimated as is done in Böttcher (2005)

$$\begin{equation} \Delta t_{\rm fl} \sim {R_m \over 8 \, \Gamma ^4 \, c} \sim 1.3 \times 10^3 \, R_{m, pc} \, \Gamma _1^{-4} \; {\rm s}. \end{equation} \tag{ 2 }$$

With this timescale, we can calculate an average flux received by the cloud as

$$\begin{eqnarray} F_m^{\rm ave} &\approx& {F_x^{\rm qu} \, d_L^2 \over c \, \Delta t_{\rm fl}} \int \limits _{R_b}^{R_m} {dx \over x^2} \approx {F_x^{\rm qu} \, d_L^2 \over c \, \Delta t_{\rm fl} \, R_b}\nonumber\\ && \sim 5.3 \times 10^{13} \; R_{m, pc}^{-1} \, \Gamma _1^4 \, R_{b, 16}^{-1} \; {\rm erg \; cm}^{-2} \; {\rm s}^{-1}, \end{eqnarray} \tag{ 3 }$$

where we have used an estimate of the quiescent synchrotron X-ray flux from 1ES 1959+650 of $F_x^{\rm qu} = 1.0 \times 10^{-10}$ erg cm⁻² s⁻¹ (similar to the average found from the Swift XRT observations), and R_b = 10¹⁶ R_{b, 16} cm is the radius of the newly emerged jet component (the "blob"). For a characteristic density of the clouds of n_c ∼ 10⁶ cm⁻³, this leads to an ionization parameter of

$$\begin{equation} \xi _{\rm ion} \equiv {F_m^{\rm ave} \over n_c} \sim 5 \times 10^7 \; {\rm erg \; cm \; s}^{-1}. \end{equation} \tag{ 4 }$$

This implies that any dust is expected to be destroyed, and that all gas is to be highly ionized by the impinging X-ray emission. If the cloud is thermally reprocessing this flux, this would require an equilibrium temperature T_equi ∼ 36, 000 K, which also requires the gas to be highly ionized. We therefore conclude that the most likely mode of reprocessing the accumulated jet synchrotron emission is Compton reflection of free electrons in the highly ionized cloud. It is not necessary that this is the only cloud within the vicinity of the blazar, maintaining the possibility that there is additional neutral absorbing gas surrounding the blazar, as found in Furniss et al. (2013). We only assume that a blob moving relativistically down the jet passes sufficiently close to a cloud that radiation from the blob can temporarily ionize the cloud, so that for a short time the blob emission is reprocessed via Compton reflection off free electrons in the cloud. This ionized gas might also act as a shield to molecular gas such as CO, as predicted in photodissociation region models (e.g., Tielenn & Hollenbach 1985; Krumholz et al. 2009; Glover & Clark 2012).

The characteristic photon frequency of jet synchrotron photons from 1ES 1959+650 is ν_sy ∼ 10¹⁷ Hz (see Figure 5), corresponding to a normalized photon energy of _sy ≡ hν_sy/(m_ec²) ∼ 10⁻². Upon Compton reflection by the cloud, this will be boosted in the jet rest frame to $\epsilon ^{\prime }_{\rm sy} \sim 0.1 \, \Gamma _1$. Therefore, any relativistic electrons (with γ_e ≳ 10) will interact with these reflected photons in the Klein–Nishina regime, resulting in strongly suppressed Compton scattering, making the production of a gamma-ray flare after the emergence of a new blob within the jet unlikely in a purely leptonic scenario.

It can be seen that the bulk Lorentz factor of the blob is a critical unknown parameter in this model, with key derived parameters strongly depending on it. The observed synchrotron peak frequency of 1ES 1959+650 is ∼1 × 10¹⁷ Hz. Therefore, even without any blueshifting from bulk motion, Compton scattering in the Thomson regime happens up to _C ∼ 1/_sy ∼ 100, corresponding to ∼50 MeV. Therefore, even for Γ = 1, a synchrotron mirror scenario would not efficiently produce VHE γ-rays via Compton scattering.

Relativistic protons with Lorentz factors of γ_p ≳ 6 × 10³, on the other hand, can interact with the reflected photons through pion production processes. Following the analysis in Böttcher (2005) and using the average flux $F_m^{\rm ave}$, we find that producing a VHE γ-ray flare with a luminosity of L_VHE ∼ 1.5 × 10⁴⁵ erg s⁻¹ requires a total number density n_p of relativistic protons

$$\begin{equation} n_p \sim 1.4 \times 10^5 \, R_{m, pc} \, \Gamma _1^{-12} \, \tau _{-1}^{-1} \, R_{b, 16}^{-2} \; {\rm cm}^{-3}, \end{equation} \tag{ 5 }$$

where τ_m = 0.1 τ₋₁ is the fraction of incident flux reflected by the cloud. We have assumed that the protons have a power-law distribution in energy $N(\gamma _p) \propto \gamma _p^{-2}$ with a low-energy cutoff at γ_{p, min} = Γ. This corresponds to a total (co-moving frame) energy in relativistic protons in the blob of

$$\begin{equation} E^{\prime }_{b, p} \sim 3.2 \times 10^{49} \, R_{m, pc} \, \Gamma _1^{-11} \, \tau _{-1}^{-1} \, R_{b, 16} \; {\rm erg}, \end{equation} \tag{ 6 }$$

and a kinetic power in the jet, carried by relativistic protons, of

$$\begin{equation} L_p \sim 7.3 \times 10^{45} \, R_{m, pc} \, \Gamma _1^{-9} \, \tau _{-1}^{-1} \, f \; {\rm erg \; s}^{-1}, \end{equation} \tag{ 7 }$$

where f is the filling factor of the jet, i.e., the fraction of the jet length occupied by plasma containing relativistic protons.

The power requirement in Equation (7) is quite moderate if one allows for a plausible filling factor f ≲ 0.1. Also, note the extremely strong dependence of the estimates in Equations (5)–(7) on the Lorentz factor. A value of Γ just slightly above 10 will reduce all energy requirements to very reasonable values, corresponding to a population of relativistic protons with energies between 10 and 100 TeV.