© ESO 2019

1. Introduction

3C 279 (redshift z  =  0.536, Burbidge & Rosenberg 1965; Marziani et al. 1996, RA_J2000  =  12^h56^m11.1^s, Dec_J2000 = −05^d47^m22^s) belongs to the class of flat spectrum radio quasars (FSRQs) that are characterized by strong variability in all energy bands from radio to γ-rays, and broad emission lines (equivalent width > 5 Å) in the optical spectrum signifying the existence of a broad-line region (BLR). FSRQs belong to the blazar class of active galactic nuclei, and their jets are closely aligned with the line of sight (Blandford & Rees 1974) resulting in strongly Doppler-boosted emission. Spectral energy distributions (SEDs) of FSRQs exhibit two broad, non-thermal components. The low-energy component peaks in the infrared and is attributed to electron synchrotron emission. In leptonic scenarios, the high-energy component, which peaks below the GeV regime, is attributed to inverse Compton (IC) emission of the same electrons scattering off ambient, soft photon fields. Such soft photon fields can be the synchrotron emission (synchrotron-self Compton, or SSC), photons from the accretion disk (IC/Disk), the broad-line region (IC/BLR), or the infrared emission of the dusty torus (IC/DT). In lepto-hadronic models, the high-energy spectral component is attributed to processes involving highly relativistic protons, such as proton synchrotron, or secondary emission from photo-meson production. The latter includes synchrotron emission from charged pions, muons, and the resulting secondary electrons and positrons. For a review of these processes see, for example, Böttcher (2007).

While FSRQs are bright in the high-energy (HE, 100 MeV <  E <  100 GeV) γ-ray domain, they are much fainter at very-high-energy (VHE, E >  100 GeV) γ-rays for a number of reasons. Firstly, the low peak energy around the lower end of the HE γ-ray domain might indicate a low maximum particle Lorentz factor, implying emission well below the VHE regime. Secondly, if the γ-rays are produced within ∼0.1 pc from the central supermassive black hole, any VHE emission would be strongly attenuated by the BLR photon field. Observations of VHE emission will therefore allow one to significantly constrain the minimum distance of the emission region from the black hole as the intrinsic absorption by the BLR cannot be too severe. Thirdly, FSRQs are found at rather large cosmological redshifts, with the closest VHE-detected FSRQ at z = 0.189 (PKS 0736+017, Cerruti et al. 2017). Hence, attenuation of VHE γ-rays by the extragalactic background light (EBL) will also reduce the detectable γ-ray flux.

3C 279 was detected at VHE γ-rays with MAGIC in 2006 (MAGIC Collaboration 2008) and 2007 (Aleksic et al. 2011) during bright optical flares. However, it has not been detected at VHE γ-rays since then (H.E.S.S. Collaboration 2014; Aleksic et al. 2014; Archambault et al. 2016). In the HE γ-ray regime, 3C 279 was detected with both EGRET (Hartman et al. 1999) and Fermi-LAT (Acero et al. 2015). Due to the ongoing monitoring of Fermi-LAT, several flares of 3C 279 have been observed in the last years, a few of which have been subject to follow-up observations with Cherenkov experiments.

In April 2014 and June 2015, 3C 279 exhibited strong outbursts in the HE γ-ray band with integrated fluxes exceeding 10⁻⁵ ph cm⁻² s⁻¹ on timescales of a few hours (Hayashida et al. 2015; Paliya 2015). Both flares were observed with Fermi-LAT in pointing mode, that is instead of the usual survey mode, the satellite was pointed towards 3C 279 to increase the exposure. In the 2015 event, this resulted in the detection of very fast variability on the order of a few minutes (Ackermann et al. 2016) on top of the longer-term (several hours) evolution of the event. Both of these events have been followed up with the High Energy Stereoscopic System (H.E.S.S.), and the results are reported here. While there is no detection in VHE γ-rays in 2014, the 2015 observation has resulted in a significant detection.

This paper is organized as follows: Sect. 2 describes the analysis of the H.E.S.S. observations of both flares. Given the H.E.S.S. detection in 2015, the analysis of a multiwavelength data set of that event is presented in Sect. 3. Sections 4 and 5 are devoted to a discussion and interpretation of both events based on various models, with an emphasis placed on the 2015 event. Limits on Lorentz invariance violations (LIV) are derived in Sect. 6. The results are summarized in Sect. 7.

Throughout the paper a Lambda cold dark matter cosmology is used with H₀ = 69.6 km s⁻¹ Mpc⁻¹, Ω_M = 0.286, and Ω_λ = 0.714 (e.g., Bennet et al. 2014). The resulting luminosity distance of 3C 279 is d_L = 3.11 Gpc.

2. H.E.S.S. data analysis

H.E.S.S. is located in the Khomas Highland in Namibia at about 1800 m above sea level. It is an array of five Imaging Atmospheric Cherenkov Telescopes, with four telescopes (CT 1−4) with 107 m² mirror area arranged in a square of 120 m side-length and one telescope (CT 5) with 614 m² mirror area in the center of the array. Observations are carried out in individual runs of typically 28 min duration. For point sources, such as 3C 279, the array observes in wobble mode, meaning with alternating offsets to the source in right ascension and declination between runs for improved background subtraction. While the array operates in stereo mode – all telescopes point at the same sky coordinate – the analysis can be performed for different array layouts depending on the demands of the observed source. A stereo analysis requires that Cherenkov emission be detected by at least two telescopes, while a mono analysis considers photons detected by CT 5. A mono analysis with CT 5 typically provides a lower energy threshold compared to analyses including CT 1−4 owing to the larger mirror area. The main analysis is performed using the MODEL analysis chain (de Naurois & Rolland 2009; Holler et al. 2015). It is cross-checked with an independent calibration chain and the analysis software IMPACT (Parsons & Hinton 2014; Parsons et al. 2015).

In 2014, H.E.S.S. observed 3C 279 with the full array over three consecutive nights between April 2 and April 4 (MJD 56749–56751). A mono analysis has been conducted with very loose cuts¹ (H.E.S.S. Collaboration 2017) resulting in an energy threshold of 66 GeV. Seven observation runs passed the quality selection (H.E.S.S. Collaboration 2006), resulting in 2.6 h of acceptance-corrected observation time, and yielding a 3.6σ significance following Li & Ma (1983). Differential upper limits (99% confidence level) have been derived following Feldman & Cousins (1998) assuming a photon index of 4. The index has been motivated by the detection spectrum of MAGIC Collaboration (2008). The upper limits are shown in Fig. 1.

Fig. 1. Observed H.E.S.S. photon spectra for six data sets as labeled. Arrows mark upper limits (99% confidence level). The gray butterfly is the 1σ statistical uncertainty band of the 2015/Night 2 data set. Error bars are statistical only. The second label gives the telescope participation and the analysis used.

Observations in 2015 were conducted in five nights between June 15 and June 21 (MJD 57188–57194) with changing array configurations. During the first night, June 15 (MJD 57188.7–57188.9, “Night 1”), CT 5 was unavailable, and a stereo analysis with loose cuts² (H.E.S.S. Collaboration 2006) has been conducted on events recorded by CT 1−4 yielding an energy threshold of 216 GeV. Quality selection has resulted in six observation runs for the analysis with 2.2 h of acceptance corrected observation time and a significance of 1.5σ. As for 2014, differential upper limits have been computed with a photon index of 4, cf. Fig. 1. Additionally, an integrated upper limit above 200 GeV has been computed, which is shown in the lightcurve in Fig. 2a.

Fig. 2. Observed multiwavelength lightcurves. (a) H.E.S.S. lightcurve derived above an energy threshold of 200 GeV in night-wise time bins with array configuration as indicated. Arrows mark upper limits (99% confidence level). (b) Fermi-LAT lightcurve integrated above 100 MeV in 3 h bins. Gray arrows mark upper limits (95% confidence level). (c) HE γ-ray photon index measured with Fermi-LAT in 3 h bins. (d) Swift-XRT lightcurve integrated between 2 and 10 keV for individual pointings. (e) Optical R band lightcurve from ATOM and SMARTS for individual pointings. (f) Spectral index between the J and B band using SMARTS observations for individual pointings. In all panels, only statistical error bars are shown.

During the second night of observations, June 16 (MJD 57189.7–57189.9, “Night 2”), CT 5 was available, and a mono analysis has been conducted with very loose cuts and an energy threshold of 66 GeV. Quality selection has led to seven observation runs for the analysis with 2.2 h of acceptance corrected observation time, resulting in a detection with 8.7σ significance. The spectrum has been modeled assuming a power-law of the form

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}E}=N_0{\left(\frac{E}{E_0}\right)}^{-\mathrm{\Gamma }}\end{array}$$(1)

with normalization N₀  =  (2.5  ±  0.2_stat  ±  0.5_sys)  ×  10⁻⁹ cm⁻² s⁻¹TeV⁻¹, photon index Γ = 4.2  ±  0.3_stat  ±  0.2_sys, and decorrelation energy E₀ = 98 GeV; see also Table 1. The systematic errors have been derived following H.E.S.S. Collaboration (2017). The spectrum is shown as the gray butterfly (1σ statistical uncertainty band), points (> 2σ significance level) and arrows (99% confidence upper limits) in Fig. 1. There is no indication for curvature as the goodness-of-fit probability of the power-law spectrum is p = 0.82. In the following, H.E.S.S. data points that have been corrected for EBL absorption using the EBL model of Franceschini et al. (2008), are used.

The average flux above an energy threshold³ of 200 GeV equals (7.6 ± 0.7_stat ± 1.5_sys) × 10⁻¹² cm⁻² s⁻¹, and is shown in Fig. 2a. A zoom into Night 2 is shown in Fig. 3a using run-wise time bins. In order to be comparable to the results of MAGIC in 2006 and 2007 (MAGIC Collaboration 2008; Aleksic et al. 2011), here the lightcurve is derived above an energy threshold of 100 GeV. The average flux is (6.5 ± 0.6_stat ± 1.3_sys) × 10⁻¹¹ cm⁻² s⁻¹, which is a factor ∼10 less than the flux during the MAGIC detection in 2006 (MAGIC Collaboration 2008). There is no indication for statistically significant variations in this lightcurve, as a constant flux has a probability of p = 0.39 (χ²/ndf = 7.6/6).

Fig. 3. Observed multiwavelength lightcurves zoomed in on Night 2. (a) H.E.S.S. lightcurve derived above an energy threshold of 100 GeV in run-wise time bins. (b) Fermi-LAT lightcurve integrated above 100 MeV in 3 h bins. (c) HE γ-ray photon index measured with Fermi-LAT in 3 h bins. (d) Swift-XRT lightcurve integrated between 2 and 10 keV for individual pointings. (e) Optical R band lightcurve from ATOM and SMARTS for individual pointings. In all panels, only statistical error bars are shown, while horizontal bars mark the observation time.

Observations on June 17 (MJD 57190.7344–57190.8569, “Night 3”) were conducted using only CT 1, 3 and 4. Six runs passed the quality selection, and a stereo analysis with loose cuts resulted in a significance of −0.6σ in 2.3 h of acceptance corrected observation time. The differential upper limit spectrum (photon index 4) is shown in Fig. 1, while the integrated upper limit above an energy threshold of 200 GeV is shown in Fig. 2a.

On June 18 (MJD 57191.7819–57191.9193, “Night 4”) all five telescopes participated in the observations. However, only two of the five conducted runs passed the CT 5 quality selection, which is why a stereo analysis with loose cuts has been done on all five runs with only the small telescopes. The analysis resulted in a significance of −2.0σ in 1.7 h of acceptance corrected observation time. The differential upper limit spectrum is shown in Fig. 1 and was computed with a photon index of 4, while the integrated upper limit above an energy threshold of 200 GeV is given in Fig. 2a.

Two more runs were taken on June 20 (MJD 57193.8339–57193.8740, “Night 5”) with all five telescopes. However, as in Night 4, the data recorded with CT 5 did not pass the quality selection. Hence again a stereo analysis with loose cuts has been performed on the data recorded with the small telescopes. Due to moon constraints the observations started relatively late, resulting in elevations of less than 52°. This explains the high energy threshold of more than 400 GeV in this night. The significance is −0.3σ in 0.7 h of acceptance corrected observation time. As before, the differential upper limit spectrum (photon index 4) is shown in Fig. 1, while the integrated upper limit above an energy threshold⁴ of 200 GeV is shown in Fig. 2a).

While the lightcurve shown in Fig. 2a may be suggestive of variability, the upper limits and the flux point have been achieved with different array configurations. An analysis of Night 2 using only the data from CT 1−4 results in no detection with an integrated upper limit comparable to the other nights. As the multiwavelength flare subsided after Night 2, and no further detections were achieved with H.E.S.S. after that night, the following discussion will focus on Nights 1 and 2 only.

3. Multiwavelength observations of the 2015 flare

In Figs. 2 and 3 lightcurves at different wavelengths of the 2015 flare are shown. The analyses are presented below.

3.1. HE γ-ray data

For the HE band, data taken with the the Large Area Telescope (Atwood et al. 2009, LAT) on-board the Fermi satellite have been analyzed. The Fermi-LAT analysis has been carried out using the Science Tool version 10.0.5 and Instrument Response Functions (IRFs) P8R2_SOURCE_V6. Data have been analyzed first on a 28 day interval, from MJD 57174 to MJD 57202 using a Binned Analysis method (Mattox et al. 1996) on a square region of interest of 30° side length and an energy range going from 100 MeV to 300 GeV. Nearby sources have been modeled using the 3FGL catalog (Acero et al. 2015) up to a radial distance from the central source of 25°. The spectral parameters of these background sources are kept free if they are within a circle of 5° from the position of 3C 279. In the annulus with angular distances between 5° and 15° only the flux normalization is left free to vary. According to the recommendations of the Fermi-LAT collaboration, the background models iso_P8R2_SOURCE_V6_v06.txt (isotropic) and gll_iem_v06.fit (galactic)⁵ are used with their normalization fit to the data.

The lightcurve and spectra for 3C 279 are obtained by fixing all the background sources in the best fit model obtained from the 28-day time interval, leaving only the spectral parameters for 3C 279 free to vary. Due to the very high level of photon counts available with Fermi-LAT for this event, it is possible firstly to perform a detailed 3 h binned lightcurve of the source near the peak of the emission shown in Figs. 2b and 3b along with the photon index in Figs. 2c and 3c, and secondly to compute the HE γ-ray spectrum in time intervals strictly simultaneous with the first and second night of the H.E.S.S. observations. In order to create a self-consistent model of the evolution of the flare (see Sect. 5.4) two more spectra are produced, namely for the “Preflare” time frame and the “Maximum” of the Fermi-LAT lightcurve between Night 1 and Night 2. The precise integration times are given in Table 1. For the calculation of the Fermi-LAT SED points, a likelihood fit has been performed in the designated energy range, with all free parameters fixed to the best power-law fit values except the normalization of 3C 279. As for lightcurves, a flux point has been computed in case the significance in the bin is above 3σ, a 95% upper limit has been calculated otherwise, assuming the best-fit power-law photon index over the entire energy range.

In the 3FGL catalogue the HE spectrum is better described by a log-parabola function of the form

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}E}=N_0{\left(\frac{E}{E_0}\right)}^{-(\mathrm{\Gamma }+\beta log\frac{E}{E_0})}\end{array}$$(2)

with the curvature parameter β. In the short time intervals of the observations considered here, only for the Maximum time frame a curved spectrum is preferred on a 4σ significance level over a power-law. The fit parameters are as follows: N₀ = 31 ± 2 × 10⁻⁵ ph cm⁻² s⁻¹ GeV⁻¹, Γ_LAT = 1.96 ± 0.05, and β_LAT = 0.12 ± 0.03 at an energy scale E₀ = 0.342 GeV. The best fit spectral values using a power-law, Eq. (1), are reported in Table 1.

3.2. X-ray data

The Neil Gehrels Swift observatory (Gehrels et al. 2004) includes three instruments: the Burst Alert Telescope (BAT, Barthelmy et al. 2005), the X-ray Telescope (XRT, Burrows et al. 2005) and the Ultraviolet/Optical Telescope (UVOT, Roming et al. 2005). These three instruments provide coverage of the following energy ranges: 5−150 keV (BAT), 0.3−10 keV (XRT), and in six optical and ultraviolet filters in the 170−600 nm wavelength range (UVOT).

XRT data collected in 2015, with Observation Ids 00035019171-00035019188, have been analyzed using version 6.21 of the HEASOFT package⁶. Data calibration has been performed using the xrtpipeline procedure and spectral fitting of each single observation has been performed with the XSPEC software (Arnaud 1996). For the fitting, all observations have been binned to have at least 30 counts per bin and each single observation has been fit with a single power-law model with a Galactic absorption value of N_H = 2.01 × 10²⁰ cm⁻² (Kalberla et al. 2005) set as a frozen parameter.

The only strictly simultaneous Swift observation was during the Maximum time frame. For Night 1 and Night 2, observations have been chosen that were conducted close to the time frames defined in Table 1. The respective Observation IDs, as well as observation times are summarized in Table 2, while the spectral results are given in Table 1. The lightcurve is shown in Fig. 2d and zoom in on Night 2 in Fig. 3d.

3.3. UV/Optical/IR data

Simultaneously with XRT, 3C 279 was monitored in the ultraviolet and optical bands with the UVOT instrument. Observations were taken in six filters: UVW2 (192.8 nm), UVM2 (224.6 nm), UVW1 (260.0 nm), U (346.5 nm), B (439.2 nm), and V (546.8 nm) (Poole et al. 2008). Magnitudes and fluxes have been calculated using uvotsource including all photons from a circular region with radius 5″. In order to determine the background, a circular region with a radius of 10″ located near the source area has been selected. All data points are corrected for dust absorption using the reddening E(B − V) = 0.0245 mag (Schlafly & Finkbeiner 2011) and the ratios of the extinction to reddening, A_λ/E(B − V) (Giommi et al. 2006). Unfortunately, only one UVOT filter was used per Swift pointing (see Table 2) during the flare. Hence, while the resulting fluxes are used in the SED in Fig. 4, no lightcurve is shown in Fig. 2.

Fig. 4. Observed multiwavelength SED for the considered time frames with black dots for the Preflare time frame, red filled squares for Night 1, green open squares for the Maximum, and blue diamonds for Night 2. The γ-ray data have been corrected for EBL absorption using the model by Franceschini et al. (2008). The solid lines show a power-law interpolation for the X-ray to γ-ray spectrum, as described in the text.

The Automatic Telescope for Optical Monitoring (ATOM, Hauser et al. 2004) is a 75 cm optical telescope located at the H.E.S.S. site in Namibia. Since 2005, it has monitored around 300 γ-ray emitters and provides optical data for H.E.S.S. observations. In 2015, 3C 279 was monitored with ATOM in the R-band from March until August. Following a rise in flux in June and coinciding with the H.E.S.S. Target-of-Opportunity observations, coverage was increased to up to 20 exposures per night, evenly spread during the time interval from 17h30 to 21h00 UTC. The flux of each observation has been derived using differential photometry using six secondary standard stars from Gonzalez-Perez et al. (2001) in the same field-of-view. The data points have been extinction-corrected similar to the UVOT data.

SMARTS (Small and Moderate Aperture Research Telescope System) is an optical and infrared telescope dedicated for observations of Fermi-LAT blazars, visible from the SMARTS site in Chile (Bonning et al. 2012). 3C 279 has been monitored with the instrument regularly since May 2008. In this paper, the observations collected for the blazar in the season of 2015 in the B, V, R, and J bands have been analyzed. SMARTS data have been corrected for extinction using the corresponding band extinctions from the Galactic Dust Reddening and Extinction Service⁷.

The R-band lightcurve is shown in Fig. 2e, while the spectral index between the B and J band, calculated as

$$\begin{array}{c}\hfill {\alpha }_{J-B}=\frac{log\nu F_J-log\nu F_B}{log{\nu }_J-log{\nu }_B},\end{array}$$(3)

is shown in Fig. 2f. Here, νF_J and νF_B are the energy fluxes in the J and B band, respectively, while ν_J and ν_B are the respective central frequencies of the filters. A zoom-in on the R-band fluxes of Night 2 is shown in Fig. 3e.

3.4. Discussion

The HE γ-ray flux, cf. Fig. 2b, increases by roughly a factor 6 from the Preflare period to Night 1, followed by another increase by a factor ∼3. The maximum is, hence, a factor ∼20 above the Preflare value. Night 2 is a factor ∼4 below the maximum and about 30% below the Night 1 flux.

The X-ray flux, cf. Fig. 2d, increases by a factor ∼2 from Night 1 to the Maximum, and drops subsequently by a factor ∼3.5. These are similar to the ratios of the HE γ-ray lightcurve and indicate a roughly simultaneous variation of the two bands.

The optical R-band flux rises by about 40% from the Preflare to Night 1, and is at a similar value in Night 2, as is shown in Fig. 2e. The detailed lightcurves from ATOM, as given in Fig. 3e, indicate minor intranight fluctuations. However, the average value is a good indicator of the optical flux state across the observation window.

Lightcurves are typically exploited to derive a characteristic timescale of a flaring event. For the 2015 flare, Ackermann et al. (2016) derived a flux doubling timescale of less than 5 min during the Maximum time frame. However, as the flare bracketed by Nights 1 and 2 lasts for roughly a day, a timescale on the order of minutes is not representative of the whole event. From the HE γ-ray lightcurve in Fig. 2b, the rise time from the low-point around Night 1 to the Maximum is about 9 h. The subsequent decay is well described by an exponential function, if the small fluctuations on top of the trend are disregarded. An exponential decay is expected from particle cooling, or if the particles leave the emission region on an energy independent timescale. Performing an exponential fit to the decaying lightcurve, one obtains a timescale of ∼9 h. Hence, this value is considered as the characteristic timescale of the event.

The observed multiwavelength SEDs are shown in Fig. 4 for the time frames defined in Table 1. In cases where multiple observations are available within a time frame, the data have been averaged. The spectral parameters of individual frequency ranges are important for modeling purposes, since they reveal information about the underlying particle distribution.

The high fluxes during the flaring event allow a precise determination of the spectral index in the HE γ-ray band in the 3 h time bins, as shown in Fig. 2c. During the flaring event the index is ∼2.2, and hardens significantly to ∼2.0 during the Maximum between Night 1 and 2 (see also Fig. 3c). Afterwards the index softens while the flux returns to the quiescence level. At this flux level, the error on the index becomes large for 3 h time bins, and no further conclusions can be drawn as the evolution of the index. The specific parameters for the averaged spectra shown in Fig. 4 are listed in Table 1.

The X-ray spectrum changes significantly during the flare, as given in Table 1. The spectrum hardens from Night 1 to the Maximum, and softens to Night 2 with the spectrum of Night 2 being even softer than the one in Night 1. Extrapolating the X-ray spectra towards the γ-ray domain would overpredict the γ-ray fluxes in all time frames.

Hence, the broad range of frequencies between the Swift-XRT and Fermi-LAT spectrum (the explicit energy ranges are given in Table 3) has been interpolated. It is assumed that the frequency range can be fit by a power-law with spectral index α, that is the energy flux is described by νF_ν ∝ ν^α with the spectral flux density F_ν. The resulting indices are reported in Table 3 and the interpolation is shown in Fig. 4. The index is positive and constant within errors during the flare with α ∼ 0.44. Unfortunately, there is no information on the Preflare time frame. The indices of the interpolation are softer than the X-ray spectral indices⁸. While the X-ray spectra themselves are compatible with simple power-laws, their spectral points and the interpolation lines in Fig. 4 are suggestive of a break above a few keV.

The indices in the optical energy range between the J and the B band, given in Table 3 and shown in Fig. 2f, are derived from the SMARTS observations as described in the previous section. The spectrum softens significantly from the Preflare time frame to the flare, but is roughly constant during Nights 1 and 2. Swift-UVOT observations during the Maximum and Night 2 time frames utilized the UVW2 filter. As can be seen in Fig. 4, their fluxes are compatible, and the Night 2 data point agrees well with an extrapolation of the other optical points. This indicates that the optical to UV flux may have been constant during the maximum of the flare. Another possibility could be that the flux in the optical band increased, but the spectrum softened in order to preserve the UV flux.

4. The flare in April 2014

The multiwavelength data of the flare in 2014 were analyzed, modeled and discussed by Paliya et al. (2015) and Hayashida et al. (2015). Paliya et al. (2015) provide a 3 h-binned HE lightcurve obtained with Fermi-LAT. This allows one to get the HE γ-ray fluxes during the H.E.S.S. observation window. They are ∼3 × 10⁻⁶ ph cm⁻² s⁻¹, ∼4 × 10⁻⁶ ph cm⁻² s⁻¹, and ∼4 × 10⁻⁶ ph cm⁻² s⁻¹, respectively. These fluxes coincide with low-points in the lightcurve between separated peaks, similar to Night 1 and Night 2 of the 2015 campaign (cf. Fig. 2). The HE fluxes in 2014 are a factor 2 to 3 lower than during Night 1 and 2 of 2015, which explains the non-detection at VHE energies.

Paliya et al. (2015) produced a HE spectrum integrated over 6 days since MJD 56749, which encompasses the H.E.S.S. observations. The average spectrum is significantly curved with photon index Γ_LAT = 2.05 ± 0.05 and curvature β_LAT = 0.13 ± 0.03⁹. These parameters are compatible with the parameters obtained in Sect. 3.1 for the Maximum time frame of 2015. The normalization for the Paliya et al. (2015) spectrum is N₀ = 5.0 × 10⁻⁶ ph cm⁻² s⁻¹ GeV⁻¹, about a factor 5 below the normalization of the Maximum time frame in 2015. Extrapolating the Paliya et al. (2015) spectrum to 100 GeV (using the corrected value for β_LAT) one obtains an energy flux of 6.7 × 10⁻¹² erg cm⁻² s⁻¹, which is below the H.E.S.S. upper limit at that energy (cf. Fig. 1).

Hayashida et al. (2015) derived a HE spectrum for a 6 h time period around the maximum flux (integration time: MJD 56750.210–56750.477), which is between the first and second night of the H.E.S.S. observations in that year. The derived HE spectrum is compatible with a power-law. The parameters are Γ_LAT = 2.16 ± 0.06, and N₀ = 1.3 × 10⁻⁵ ph cm⁻² s⁻¹ GeV⁻¹, which are similar to the parameters obtained for Night 2 in 2015. Hence, a detection at VHE may have been possible during the peak flux in 2014.

Paliya et al. (2015) and Hayashida et al. (2015) used leptonic one-zone models using different combinations of SSC, IC/BLR and IC/DT emission for the high-energy peak. The H.E.S.S. upper limits cannot constrain the models.

5. The flare in June 2015

The significant detection of the 3C 279 flare with H.E.S.S. in 2015 gives important constraints on the parameter space. These constraints are discussed below, and time-dependent leptonic and lepto-hadronic one-zone models are tested to account for the variability. Most notably, the combined fit of the Fermi-LAT and H.E.S.S. spectra in Night 2 provides strong constraints on the absorption of γ-rays, which can be used to constrain the minimum distance of the emission region to the black hole. This is presented first, followed by a brief description of the prevalent thermal photon fields surrounding the jet, which will be used for both modeling attempts.

5.1. Minimum distance of the emission region from the black hole

The contemporaneous data of Fermi-LAT and H.E.S.S. enable the search for absorption features caused by pair production of γ-rays with photons of the BLR. The latter is derived following the model of Finke (2016), which is motivated by reverberation mapping and assumes that accretion disc radiation is absorbed by the BLR clouds and re-emitted as monochromatic lines at fixed distances from the black hole. The approach here closely follows the method introduced by Meyer et al. (2019), who used Fermi-LAT data of six bright FSRQs to search for absorption features.

Two geometries of the BLR are implemented in the model. In the shell geometry, BLR photons are emitted in infinitesimally thin shells around the black hole, whereas in the ring geometry, the BLR photons originate from thin rings orthogonal to the jet axis. The model includes emission lines from Lyϵ to Hα but neglects any contribution from the thermal continuum. Each line has an associated luminosity and is emitted in a shell or a ring at a fixed distance (see Table 5 in Finke 2016). As input the model requires the black hole mass, M_•, and the luminosity of the Hβ line, L(Hβ). For 3C 279, log₁₀(M_•/M_⊙) = 8.28 with the solar mass M_⊙, and L(Hβ) = 1.7 × 10⁴³ erg s⁻¹ are adopted (Liu et al. 2006). Using the relations summarized in Finke (2016) between L(Hβ) and L(5100 Å), as well as between L(5100 Å) and the radius of the Hβ emitting shell together with Table 5 of Finke (2016), the radius of the Lyα emitting shell, R_Lyα ∼ 7.6 × 10¹⁶ cm, is obtained. The Lyα luminosity is the highest in the model (a factor of 12 higher than L(Hβ)) and is therefore responsible for most of the absorption. The values for R_Lyα and the Lyα luminosity are broadly consistent with typical values obtained from reverberation mapping (Kaspi et al. 2007; Bentz et al. 2009; Meyer et al. 2019) The resulting optical depths, τ_γγ(r, E), for both geometries and different distances r of the emission region from the central black hole are shown as a function of the γ-ray energy in Fig. 5. The shell geometry generally results in higher values of the optical depth (compare also Fig. 14 in Finke 2016). Nevertheless, the optical depths are still lower compared to predictions of more sophisticated BLR models that include continuum emission (e.g., Abolmasov & Poutanen 2017, see also the discussion in Meyer et al. 2019). In that sense, constraints on the minimum distance between the γ-ray emitting region and the central black hole can be regarded as conservative.

The distance r is constrained by simultaneously fitting the Fermi-LAT and H.E.S.S. data, both corrected for the EBL influence following Franceschini et al. (2008), with an intrinsic spectrum F(E) which is modified by the absorption exp(−τ_γγ) (Fig. 6, left). The EBL model of Franceschini et al. (2008) is in good agreement with other EBL models and with lower limits derived from galaxy number counts (see Dwek & Krennrich 2013, for a review). Since a spectral cut-off due to absorption is degenerate with a cut-off of the intrinsic spectrum, different intrinsic spectral shapes, namely a log-parabola (LP), a power law with sub-exponential cut-off (SEPL) and a broken power law (BPL) are tested. For each combination of intrinsic spectrum and assumed BLR geometry (ring or shell), the parameters of the intrinsic spectrum and r are optimized. This is done using a maximum likelihood optimization, where the likelihood of each Fermi-LAT and H.E.S.S. spectral flux point is approximated with a Gaussian centered on the measured flux and with a width equal to the flux uncertainty in each bin. One-sided Gaussian distributions are used in case of flux upper limits.

Fig. 5. Optical depths for γ-rays emitted along the jet axis at different distances r interacting with photons of the BLR emission lines. Left: BLR modeled with the shell geometry. The crossing of lines at low energies is due to numerical inaccuracies. Right: BLR modeled with the ring geometry. The structure in the optical depth are caused by the contributions of different emission lines to the overall optical depth.

Fig. 6. Left: best-fit spectra for the BLR ring geometry to the combined Fermi-LAT and H.E.S.S. data. Both data sets are corrected for EBL absorption following Franceschini et al. (2008). The spectral shapes do not change significantly, if the shell geometry is assumed instead. Right: likelihood profile as a function of the distance r for the different assumed intrinsic spectra and BLR geometries.

The resulting best-fit spectra for the ring geometry are shown in the left panel of Fig. 6. The best-fit values for r are around ∼11 R_Lyα for the ring geometry and around ∼10 R_Lyα in the shell geometry regardless of tested spectral shapes. The figure includes the χ² values per degrees of freedom (d.o.f.). The reduced χ² values are all above unity and the fit qualities, measured by the p-value of the χ² distribution with corresponding d.o.f., are 0.11, 0.12, 0.06 (0.01, 0.01, 0.003) for the LP, BPL, SEPL intrinsic spectra and the ring (shell) geometry, respectively. For the LP case, the dotted line additionally shows the case when r is fixed to 2 R_Lyα. For such small values of r, the BLR absorption leads to a sharp cut-off of the observed spectrum. We note that for the SEPL case, a sub-exponential cut-off is preferred by the data. A standard exponential cut-off could reproduce the Fermi-LAT data and the first two flux points obtained with H.E.S.S. but would under-predict the flux in the highest energy bin by an order of magnitude.

The right panel of Fig. 6 shows the profile likelihood of the fit as a function of r. It is evident from the figure that none of the fits significantly prefers the presence of an absorption feature at these large distances over the no-absorption case (which corresponds to the maximum tested distance, r ∼ 30 R_Lyα). Therefore, the maximum likelihood approach is used to derive 95% confidence lower limits on r. The lower limits are found by decreasing r until the likelihood increases by Δ ln ℒ = 2.71/2. All assumed intrinsic spectra result in roughly the same value of the limit of r ≳ 5.4 R_Lyα = 4.1 × 10¹⁷ cm. Since the optical depth is smaller for the ring geometry, the lower limit in this case relaxes to r ≳ 2.6 R_Lyα = 2.0 × 10¹⁷ cm for the LP and SEPL intrinsic spectra. The lower limit is slightly lower for the BPL spectrum, r ≳ 2.2 R_Lyα = 1.7 × 10¹⁷ cm. It should be noted that if only the Fermi-LAT data points are fit with a power law, which is then extrapolated to higher energies including BLR absorption, the flux for all HESS data points is severely under-predicted for r ≲ 7 × 10¹⁶ cm. This model does not provide a satisfactory fit to the H.E.S.S. data and is especially in tension with the highest energy H.E.S.S. data point, which it under-predicts by more than an order of magnitude. In conclusion, the emission zone is confidently placed beyond r ∼ 1.7 × 10¹⁷ cm (or 3 × 10³ Schwarzschild radii), outside the BLR.

5.2. The external photon fields

In this section, the photon fields external to the jet of 3C 279 are described. The parameters are listed in Table 4 and are used for the leptonic and lepto-hadronic models described in the next sections.

The accretion disk is modeled as a Shakura-Sunyaev disk (Shakura & Sunyaev 1973) with a luminosity L_acc = 3.0 × 10⁴⁵ erg s⁻¹, which is the average of values given in the literature (e.g., Hayashida et al. 2015; Paliya et al. 2015). The accretion disk luminosity is about 8% of the Eddington power L_edd = 3.78 × 10⁴⁶ erg s⁻¹ of black hole with mass M_bh ∼ 3 × 10⁸ M_⊙ (and references therein Hayashida et al. 2015). The inner radius of the disk is set to the innermost stable orbit of a Schwarzschild black hole, namely R_acc, in = 6 ×  R_g with the gravitational radius of the black hole R_g. The outer radius can be estimated following Netzer (2015), and marks the point where the self-gravity of the disk surpasses the gravity of the black hole leading to disk fragmentation further out. For 3C 279 this corresponds to R_acc, out ∼ 430 × R_g.

Unlike the lines, the thermal BLR parameters are not well known for 3C 279. Using the numbers from the previous section, the radius of the BLR is r_BLR = R_Lyα, and the luminosity is assumed as L_BLR = 2.3 × 10⁴⁴ erg s⁻¹. This corresponds to about 8% of the accretion disk luminosity. The given BLR luminosity contains the sum of the line luminosities plus a thermal contribution. The BLR temperature is set to T_BLR = 1.0 × 10⁴ K. Note that for the inverse Compton process the BLR line emission can be well approximated by a thermal continuum.

As the discussion in Sect. 5.1 indicates that the emission region is located beyond the BLR, its emission may be an inefficient target for the IC process. Whether the strong accretion disk radiation is a useful target field despite being strongly de-boosted, cannot be stated a priori. Therefore, we also invoke the thermal field of a dusty torus, despite the fact that there is no evidence of its presence in 3C 279. Using estimates from Hayashida et al. (2012), the radius of the DT becomes r_DT = 4.23 × 10¹⁸ cm, while the luminosity in this case is assumed to be 10% of the accretion disk. The temperature is assumed to be T_DT = 500 K.

5.3. Leptonic one-zone model

The leptonic one-zone model is still the standard model for blazar physics, either in the most fundamental version with synchrotron-self Compton (SSC) or in the slightly extended version with external fields, such as the accretion disk, the broad-line reagion (BLR) and the dusty torus (DT). Its advantage is the relatively low number of parameters, of which a lot can be constrained. From now on parameters marked with a prime are considered in the host galaxy frame, quantities with an asterisk are in the observer’s frame, and unmarked quantities are either in the comoving jet frame or invariant.

The parameters used for the modeling are listed in Tables 4 and 5. Proper explanations of the parameters and the description of the code are given in Appendix A. Some of the parameters can be analytically constrained, which is also described in Appendix A.

The modeling aims to reproduce the flare at the time around the H.E.S.S. observations. Hence, first the Preflare SED is reproduced with the parameters listed above the horizontal line in Table 5, followed by Night 1. Then the Maximum is modeled, after which the evolution is followed to Night 2. The variability is modeled with the following parameter changes:

$$\begin{array}{cc}& B(t)=B+\mathrm{\Delta }{B}_1(\mathrm{H}\left[t,t_{\mathrm{s}}^{\ast },t_1^{\ast }\right]+\mathrm{H}\left[t,t_{\mathrm{m}}^{\ast },t_2^{\ast }\right])\hfill \\ \hfill & +\mathrm{\Delta }{B}_2\mathrm{H}\left[t,t_1^{\ast },t_m^{\ast }\right]\hfill \end{array}$$(4)

$$\begin{array}{cc}& L_{\mathrm{inj}}^{\mathrm{e}}(t)=L_{\mathrm{inj}}^{\mathrm{e}}+\mathrm{\Delta }{L}_{\mathrm{inj},1}^{\mathrm{e}}(\mathrm{H}\left[t,t_{\mathrm{s}}^{\ast },t_1^{\ast }\right]+\mathrm{H}\left[t,t_{\mathrm{m}}^{\ast },t_2^{\ast }\right])\hfill \\ \hfill & +\mathrm{\Delta }{L}_{\mathrm{inj},2}^{\mathrm{e}}\mathrm{H}\left[t,t_1^{\ast },t_m^{\ast }\right]\hfill \end{array}$$(5)

$$\begin{array}{cc}& {\gamma }_{\min }^{\mathrm{e}}(t)={\gamma }_{\min }^{\mathrm{e}}+\mathrm{\Delta }{\gamma }_{\min }^{\mathrm{e}}\mathrm{H}\left[t,t_{\mathrm{s}}^{\ast },t_2^{\ast }\right]\hfill \end{array}$$(6)

$$\begin{array}{cc}& s^{\mathrm{e}}(t)=s^{\mathrm{e}}+\mathrm{\Delta }{s}^{\mathrm{e}}\mathrm{H}\left[t,t_{\mathrm{s}}^{\ast },t_2^{\ast }\right],\hfill \end{array}$$(7)

where $t_{\mathrm{s}}^{*}$ = MJD 57186.875 marks the beginning of the flaring event, $t_1^{\ast }$ = MJD 57188.875 marks Night 1, $t_{\mathrm{m}}^{*}$ = MJD 57189.125 the Maximum and $t_2^{\ast }$ = MJD 57189.875 Night 2. While these time steps are defined in the observer’s frame, they are properly transformed to the comoving frame in the code. The step function H[x,a,b] is 1 for a ≤ x ≤ b and 0 otherwise. Hence, the variability is modeled by 1 or 2 step-function-like changes in the parameters. The variability parameters are listed in Table 5 below the horizontal line. A reasoning for the adopted parameter changes is provided in Appendix A.

The resulting model SEDs are shown in Fig. 7. The optical regime is dominated by synchrotron photons, while the X-ray regime is mostly SSC, and the γ-ray regime is dominated by the IC/BLR process. The SEDs are reproduced well for the Preflare, Night 1 and Maximum time frames except in the X-ray domain. However, these time frames can be directly influenced by the changes in the parameters. Subsequently, the injection is returned to Night 1 levels, so the continuing evolution is given by the cooling and escape of the particles. As Night 2 is not reproduced well in the X-ray and γ-ray energy bands, the chosen parameter set is not adequate to reproduce the decay from the Maximum to Night 2.

Fig. 7. Multiwavelength spectra and models for the four time frames: Preflare (black dots), Night 1 (red filled squares), Maximum (green open squares), and Night 2 (blue diamonds). The γ-ray data have been corrected for EBL absorption using the model by Franceschini et al. (2008). The thick solid lines mark the leptonic models. The thin lines mark spectral components for the Preflare period as labeled.

In order to improve the fit in the X-ray domain, a higher SSC flux is required. This could be achieved by a larger number of particles, which would however also increase the synchrotron and IC/BLR fluxes. This could be alleviated by reducing the magnetic field and the luminosity of the BLR. However, the latter is already close to the allowed flux from the line measurements. Increasing the magnetic field, which would in turn increase the SSC flux, would require a brighter BLR in order to preserve the Compton dominance. Additionally, this would require less particles in the emission region. As the SSC flux depends linearly on the magnetic field but quadratically on the particle density, the SSC flux would actually drop.

The bad fit of Night 2 is driven by the slow particle escape due to the large emission region. Instead of leaving the source, the particles are shifted to lower energies. This has no consequences for the optical domain, which is dominated by particles that cool quickly, explaining the good fit. However, the X-ray and HE γ-ray domains are dominated by the inverse Compton radiation of less energetic particles. In this energy regime particles have piled up, as the original ones have not yet cooled away, while further particles have reached this energy by cooling down from higher energies.

This could be alleviated by a faster escape of the particles from the emission region. As the escape is controlled by both the size of the emission region R, and the escape timescale parameter η_esc, either of them could be reduced to accelerate the escape. However, η_esc is already set to only 5.0, implying that particles remain within the emission region for only five light crossing timescales. This is already a very fast escape, as one expects some diffusion within the emission region due to the magnetic field.

Hence, reducing R is used to accelerate the escape of particles, as the constraint from the characteristic variability timescale only provides an upper limit. However, reducing R enhances the energy densities of particles and photon fields within the emission region. While this can be accommodated easily for the synchrotron and external-Compton component by reducing the number of injected particles, the SSC flux would drop, as outlined above, and therefore make the fit even worse.

Another possibility is to (additionally) increase the Doppler factor δ. As this value has a direct impact on the internal energy densities of the emission region, the parameters have to be changed considerably. However, also in this case a perfect fit is not possible under the given constraints, which is shown in Fig. A.1.

It should be noted that despite the mentioned problems, the H.E.S.S. spectrum is fit well. If the escape problem could be solved, the fit would actually be perfect in the H.E.S.S. domain as the Fermi-LAT spectra of Night 1 and 2 are similar, and so would be the models.

As mentioned above, a higher SSC flux could be achieved with a larger number of particles in the jet, while reducing the magnetic field and the external field. While reducing the BLR luminosity is not possible, the emission region could be moved to an even further distance from the black hole, where the DT photon field dominates the external contribution. In fact, parameters can be found that allow for a good fit in large parts of the spectrum, but not perfectly at all energies, cf. Fig. A.2. The main issue is again the escape of particles, but the delicate interplay of the parameters does not allow to reduce the size of the emission region in this case.

Hence, despite being able to fit the Preflare, Night 1 and Maximum time frames rather well in some cases, the subsequent decay poses a severe problem for the leptonic model. The interplay of the parameters is delicate and requires incredible fine tuning, which could not be achieved for all the details of the spectrum.

Nonetheless, it is interesting to study the resulting power output of the model shown in Fig. 7. Table 6 lists the Poynting power, proton power, electron power, and radiative power. The proton power is calculated assuming one cold proton per electron. The powers have been derived under the assumption that the bulk Lorentz factor is given by Γ_j = δ. Compared to the Eddington power of 3C 279’s black hole, L_edd = 3.78 × 10⁴⁶ erg s⁻¹, the total power is below the Eddinton limit during the Preflare and Night 1 time frames. The Maximum, and Night 2 exceed the Eddington power. By how much depends on the actual value of the mass of the black hole, which has an uncertainty of more than a factor 2 (e.g., Hayashida et al. 2015). The power output of the jet is dominated by particles and radiation, while the Poynting power is comparable to the other constituents only during the Preflare period. The total power of the jet of Night 2 could be reduced to below the Eddington limit if the emission region contains 90% pairs. Since the radiative output of the jet is already above the Eddington luminosity for the Maximum (keeping the uncertainty in M_bh in mind), even a high pair content would not be able to reduce the jet power below that threshold. It should also be noted that the model with a larger Doppler and bulk Lorentz factor (shown in Fig. A.1) results in super-Eddington jet powers, however with a smaller margin, and a high fraction of pairs may push the total jet power below the Eddington limit in this case.

5.4. Lepto-hadronic one-zone model

To go beyond the simple one-zone leptonic model, the possibility of a one-zone lepto-hadronic model is explored, following a similar set of source assumptions as made for the leptonic model. Typically, lepto-hadronic models have difficulties in reproducing fast flares owing to the long cooling timescales of protons. However, it was noted by Petropoulou et al. (2017) that small scale regions with kG magnetic fields could account for the minute-scale variability even in lepto-hadronic models. While the minute-scale variability is not a concern in the present study, it shows the principle possibility to use lepto-hadronic models to account for the 2015 flare of 3C 279.

The parameters reproducing the Preflare period are listed in Tables 4 and 7 above the horizontal line. They are explained along with a discussion of the constraints and the details of the code in Appendix B.

Again, a self-consistent reproduction of the 3C 279 spectra is attempted by changing input parameters as follows:

$$\begin{array}{cc}& L_{\mathrm{inj}}^{\mathrm{p}}(t)=L_{\mathrm{inj}}^{\mathrm{p}}+\mathrm{\Delta }{L}_{\mathrm{inj},1}^{\mathrm{p}}\mathrm{H}\left[t,t_{\mathrm{s}}^{\ast },t_1^{\ast }\right]+\mathrm{\Delta }{L}_{\mathrm{inj},2}^{\mathrm{p}}\mathrm{H}\left[t,t_1^{\ast },t_1^{\ast }\right]\hfill \end{array}$$(8)

$$\begin{array}{cc}& {\gamma }_{\max }^{\mathrm{p}}(t)={\gamma }_{\max }^{\mathrm{p}}+\mathrm{\Delta }{\gamma }_{\max }^{\mathrm{p}}\mathrm{H}\left[t,t_{\mathrm{s}}^{\ast },t_2^{\ast }\right]\hfill \end{array}$$(9)

$$\begin{array}{cc}& L_{\mathrm{inj}}^{\mathrm{e}}(t)=L_{\mathrm{inj}}^{\mathrm{e}}+\mathrm{\Delta }{L}_{\mathrm{inj}}^{\mathrm{e}}\mathrm{H}\left[t,t_{\mathrm{s}}^{\ast },t_2^{\ast }\right]\hfill \end{array}$$(10)

$$\begin{array}{cc}& s^{\mathrm{e}}(t)=s^{\mathrm{e}}+\mathrm{\Delta }{s}^{\mathrm{e}}\mathrm{H}\left[t,t_{\mathrm{s}}^{\ast },t_2^{\ast }\right],\hfill \end{array}$$(11)

where the time steps are the same as in the leptonic case. The maximum proton Lorentz factor, the electron injection luminosity, and the electron spectral index are only varied once during the flare and remain at their levels until the end of the flare. The proton injection, however, is varied until the beginning of the Maximum, with a single injection on top of that, after which it returns to Preflare levels. The variability parameters are listed in Table 7 below the horizontal line. The magnetic field is not varied, as there is no constraint on it in this case.

The four derived spectra for the lepto-hadronic model are shown in Fig. 8. The optical component is well reproduced by electron synchrotron emission in all cases. The X-ray and HE γ-rays are dominated by proton synchrotron emission, while the VHE γ-rays are influenced by muon synchrotron emission with some contributions from synchrotron emission of secondary electrons. The Preflare HE γ-ray spectrum is well matched except for the highest energy bins. Night 1 is well reproduced for both X-rays and HE γ-rays. The Maximum is well reproduced in the HE γ-rays, and X-rays. Night 2 is significantly overpredicted in both the X-rays and HE γ-rays, while the VHE γ-rays are too low.

Fig. 8. Multiwavelength spectra and models for the four considered time frames: Preflare (black dots), Night 1 (red filled squares), Maximum (green open squares), and Night 2 (blue diamonds). The γ-ray data have been corrected for EBL absorption using the model by Franceschini et al. (2008). The thick solid lines mark the lepto-hadronic photon models, while the thick dashed lines mark the neutrino models. The thin lines mark the spectral components for the Preflare period as labeled.

The main problem is, again, the slow escape of particles, which in this case is mainly the proton escape, since the electrons cool very efficiently in the strong magnetic field. However, the protons barely cool nor escape, which is why the hadronic spectral components of Night 2 even slightly exceed those of the Maximum. As before, reducing the size of the emission region R would increase particle escape, and hence reduce the flux. The increase in internal energy densities would increase the production of pions, muons and secondary electron/positron pairs. This could produce a flux that is compatible with the VHE spectrum. However, tests have revealed that the interplay of escape and cooling – while weak – has an observable effect, which makes a fit in either the X-rays or the VHE γ-rays problematic.

As discussed below, the jet power significantly exceeds the Eddington power of the black hole. This could be mitigated by increasing the bulk Lorentz and Doppler factor, as the same radiative output requires less power in the particles. A realization is shown in Fig. B.1. However, while the total jet power decreases slightly, it still significantly surpasses the Eddington power. One should also note that despite a significantly smaller source size, and a much faster escape it is still not possible to fit the X-ray and VHE γ-ray spectrum of Night 2 simultaneously, as the latter is underpredicted by the model. As the emission region is placed at the minimum distance allowed by the result of Sect. 5.1, the external fields cannot be enhanced further to allow for a larger number of muons to be produced, as their synchrotron emission is responsible for the VHE γ-ray output within this model.

It should be noted that the spectral characteristics in the X-ray domain require a rather large minimum proton Lorentz factor. This is difficult to explain through conventional acceleration processes, which expect a minimum proton Lorentz factor of ∼1.

Given the larger number of free parameters, the lepto-hadronic model is less constrained than the leptonic model. Nonetheless, the one-zone lepto-hadronic model is also not able to self-consistently reproduce the observed characteristics of the flare.

The Poynting, proton, electron/positron and radiative powers are given in Table 8. The power output in this case is dominated by the Poynting flux and the proton power, while the radiative and electron/positron powers are subdominant. The electron/positron power increases throughout the flare and even during Night 2. This comes from the ongoing injection of secondary electron/positron pairs from the muon decay, which inject highly energetic pairs that carry a large amount of power. In all cases the total power significantly exceeds the Eddington power, L_edd = 3.78 × 10⁴⁶ erg s⁻¹. The general picture does not change much for a larger bulk Lorentz factor.

The decay of pions and muons releases neutrinos, and the model neutrino spectra arriving at Earth are shown in Fig. 8 for the four specific time steps. Using IceCube’s effective area (IceCube Collaboration 2013), the detectable neutrino rates for IceCube can be calculated, and hence the potential number of detectable neutrinos from the event. Concentrating on the ∼27 h time window bounded by the H.E.S.S. observations and which covers the peak flux in the HE band, the number of detectable neutrinos is 5 × 10⁻⁴. Even if the emission region would be located within the BLR, which would result in a larger pion and muon production rate and hence a larger number of neutrinos, the rate would not increase enough in order to reach unity (Zacharias et al. 2019). Hence, no neutrino is expected to be detected by IceCube during the Maximum of the flare, and this approach cannot be used to distinguish between the leptonic and lepto-hadronic models.

5.5. Discussion

The most important result is the lower limit of the distance of the emission region from the black hole, placing it outside the BLR. This directly implies that the observed minute-scale variability in the HE γ-ray band (Ackermann et al. 2016) is not caused by an emission region encompassing the entire width of the jet. It rather points towards small emission regions or turbulent cells within a larger active region (e.g., Giannios 2013; Marscher 2014). Furthermore, it adds to the growing evidence (e.g., Zacharias et al. 2017; Costamente et al. 2018) that jets can produce γ-ray emission on large distances from the black hole.

None of the one-zone models can fully reproduce the observed characteristics of the 2015 flare in 3C 279, and the jet powers are a severe constraint for the models. The leptonic model is mostly below the Eddington power of the black hole. However, it surpasses the Eddington power during the Maximum time frame. Interestingly, during this time frame the radiative power emitted by 3C 279 is already very close to or even surpasses the Eddington limit depending on the actual mass of the black hole. This underlines the extreme nature of this flare. In the lepto-hadronic model, the Eddington power is surpassed during every time frame by a large factor. This is a common problem of proton-synchrotron models (e.g., Zdziarski & Böttcher 2015). This might be possible for a short flare, as the one described here, but it is unlikely for a longer period, such as the Preflare time frame, which resembles the ground state, where the total power of the model is ∼3 × 10⁴⁸ erg s⁻¹. The assumption that the jet power is provided dominantly by the accretion power, implies a radiative efficiency of the accretion disk of less than 10⁻³, using the bolometric disk luminosity given Table 4. This is much less than the typical radiative efficiency of accretion disks of 0.1–0.2 in active galaxies, and unlikely on long timescales (see Zdziarski & Böttcher 2015, for a detailed discussion). Hence, while the flare itself might be hadronically induced, the quiescent state is probably not.

In this work, it is assumed that the emission region is a standing shock, like a recollimation shock, within the jet and does not change its position during the flare. The bulk flow is, thus, provided by the jet material crossing the shock. However, a moving shock would cover a distance of ∼1 pc during the flare. In such a scenario, the external fields would change with time, which could explain the reduction in γ-ray flux at the end of the flare without a faster escape of particles. Additionally, more sophisticated models, such as a spine-in-sheath model (Ghisellini et al. 2005), a jets-in-jet model (Giannios 2013), a moving mirror model (Vittorini et al. 2017), and others, might provide an improved description of the observations. However, testing these possibilities is beyond the scope of this paper.

The failure of the simple leptonic one-zone model in 3C 279 has been noted before (e.g., Böttcher et al. 2009; Bottacini et al. 2016). The 2015 flare has been explicitly modeled by Bottacini et al. (2016) using stationary leptonic and hadronic one-zone models. The discussion in Bottacini et al. (2016) focuses on INTEGRAL observations conducted for ∼14 h around the peak of the γ-ray outburst. All data were integrated over this time bin, including the significant variability in all bands. Bottacini et al. (2016) conclude that their leptonic model would not produce VHE γ-ray emission. On the other hand, their hadronic model would allow for VHE γ-ray emission. This strong statement cannot be confirmed here, as the time-dependent leptonic and lepto-hadronic models allow for VHE γ-ray emission, even though a self-consistent fit could not be achieved.

6. Limits on Lorentz invariance violation

Several models of Quantum Gravity (QG) predict a violation of Lorentz invariance (LIV in the following for Lorentz invariance violation) in the form of a modified dispersion relation for photons in vacuum (see Amelino-Camelia 2013, for a general review about QG phenomenology, including modified dispersion relations). This effect should be dominant at energy scales of the order of the Planck scale (∼10¹⁹ GeV) but it is believed that it could be observed at lower energies, though linearly or quadratically suppressed. The modified dispersion relation leads to an energy-dependent speed of light. High-energy photons propagating in vacuum may either be slower (the sub-luminal case) or faster (the super-luminal case) than low energy photons. In addition, the longer the propagation distance is, the larger the time delay between photons of different energies should be.

Variable or transient astrophysical sources are then very good candidates to constrain $E_{\text{QG}}^{\text{l}}$ and $E_{\text{QG}}^{\text{q}}$, respectively the energy scales for linear and quadratic LIV effects. Stringent limits have been obtained with flaring AGN (H.E.S.S. Collaboration 2011) and pulsars (MAGIC Collaboration 2017). The best limits available so far (Vasileiou et al. 2013), obtained with Gamma-Ray Bursts GRB 090510 and GRB 090926A, are above the Planck scale for $E_{\text{QG}}^{\text{l}}$ but still need to be confirmed with more objects. 3C 279 is a prime candidate to perform such a study, owing to its high redshift and high number of VHE γ-ray photons detected during the 2015 flare (Night 2).

The search for energy-dependent time-delays was performed with the likelihood method first introduced by Martínez & Errando (2009), and adapted to take the background contribution into account as detailed by H.E.S.S. Collaboration (2015) for the case of PG 1553+113 flare of 2012. The likelihood method “compares” the arrival times of photons at high energies (here in the 300 GeV – 2 TeV range, corresponding to 185 events) with a parameterized template lightcurve obtained at low energies (100 GeV – 150 GeV, 243 events). As in H.E.S.S. Collaboration (2011, 2015), a detailed study was performed using simulations to evaluate statistical and systematic errors. Statistical errors were calibrated by generating 300 realizations of the lightcurve mimicking the data, with no lag. As in the case of PG 1553+113 flare of 2012 which shows similar characteristics of variability and background contamination, systematic errors were found to be mainly due to the low energy template lightcurve parameterization and to energy selections. Overall, no significant lag was measured and one-sided 95% confidence level limits on $E_{\text{QG}}^{\text{l}}$ and $E_{\text{QG}}^{\text{q}}$ were computed. These results are given in Table 9 for the sub- and super-luminal cases. Sub-luminal limits are also shown in Fig. 9 together with other results published by the H.E.S.S. Collaboration (H.E.S.S. Collaboration 2011, 2015, 2019). Even if less constraining than those obtained from other studies, the 3C 279 flare results will be valuable for future population studies due to the high redshift of this source.

Fig. 9. Lower limits on the sub-luminal linear (left) and quadratic (right) terms of the modified dispersion relation obtained with H.E.S.S. for several AGN as a function of redshift. 1: PKS 2155−304 (H.E.S.S. Collaboration 2011), 2: PG 1553+113 (H.E.S.S. Collaboration 2015), 3: Mrk 501 (H.E.S.S. Collaboration 2019), 4: 3C 279 (this work).

7. Summary and conclusions

The FSRQ 3C 279 underwent two major HE γ-ray outbursts in April 2014 and June 2015. These were among the brightest flares detected with Fermi-LAT, and during the 2015 flare a HE γ-ray variability timescale of ∼5 min was detected (Ackermann et al. 2016). Both flares were followed up with the VHE γ-ray experiment H.E.S.S. The observations in 2014 have not yielded a detection, and the upper limits are not particularly constraining for modeling attempts of that event. However, a significant detection of 3C 279 at VHE γ-ray has been achieved during the 2015 event with 8.7σ above an energy threshold of 66 GeV. This allows one to derive strong constraints on source parameters.

Most importantly, the VHE γ-ray spectrum along with a simultaneous HE spectrum can be used to derive the amount of absorption of the emitted γ-rays through the BLR photon field. This can be translated into an estimate of the distance of the emission region from the black hole. An elaborate analysis (Meyer et al. 2019) results in a lower limit (95% confidence level) on the distance at r ≳ 1.7 × 10¹⁷ cm, placing the emission region outside of the BLR.

In this work, using both a time-dependent leptonic, and a time-dependent lepto-hadronic one-zone model, a reproduction of the contemporaneous spectra up to VHE γ-ray energies and the hour-scale variability has been attempted. The leptonic model reproduces the data well in the optical and γ-ray bands for most time frames except for Night 2, where the spectrum is overproduced in the γ-ray component. The X-ray spectrum cannot be adequately fit. Accelerating the escape of particles after the Maximum by invoking a smaller emission region along with other parameter changes, does not improve the fit significantly.

The lepto-hadronic model faces similar difficulties, as the decrease from the Maximum to Night 2 in the high-energy component is also not well covered. Changing the parameters is also unable to provide a satisfactory fit. The number of model neutrinos is too low to be detectable by IceCube, and can therefore not be used as a discriminator of the models.

In summary, simple one-zone models cannot fully reproduce the observed characteristics of the 2015 flare in 3C 279 within the given constraints. More elaborate models are required, which are beyond the scope of this paper.

The lower limits on LIV linear and quadratic energy scales obtained in this study are comparable to those derived from other flaring AGNs observed with similar characteristics of variability and background level. The data described here will be included in the overall combination of LIV study results which is currently being prepared by the three major IACT experiments (see Nogues et al. 2017, for a preliminar study from simulated data). Due to its high redshift, 3C 279 will also be an interesting target for population studies with the future Cherenkov Telescope Array, which is expected to greatly improve the current limits on QG energy scale.

Acknowledgments

The support of the Namibian authorities and of the University of Namibia in facilitating the construction and operation of H.E.S.S. is gratefully acknowledged, as is the support by the German Ministry for Education and Research (BMBF), the Max Planck Society, the German Research Foundation (DFG), the Helmholtz Association, the Alexander von Humboldt Foundation, the French Ministry of Higher Education, Research and Innovation, the Centre National de la Recherche Scientifique (CNRS/IN2P3 and CNRS/INSU), the Commissariat á l’énergie atomique et aux énergies alternatives (CEA), the U.K. Science and Technology Facilities Council (STFC), the Knut and Alice Wallenberg Foundation, the National Science Centre, Poland grant no. 2016/22/M/ST9/00382, the South African Department of Science and Technology and National Research Foundation, the University of Namibia, the National Commission on Research, Science & Technology of Namibia (NCRST), the Austrian Federal Ministry of Education, Science and Research and the Austrian Science Fund (FWF), the Australian Research Council (ARC), the Japan Society for the Promotion of Science and by the University of Amsterdam. We appreciate the excellent work of the technical support staff in Berlin, Zeuthen, Heidelberg, Palaiseau, Paris, Saclay, Tübingen and in Namibia in the construction and operation of the equipment. This work benefited from services provided by the H.E.S.S. Virtual Organisation, supported by the national resource providers of the EGI Federation. A.W. is supported by Polish National Agency for Academic Exchange (NAWA). This paper has made use of up-to-date SMARTS optical/near-infrared lightcurves that are available at www.astro.yale.edu/smarts/glast/home.php. This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

References

Appendix A: Leptonic code and model constraints

The time-dependent leptonic code used in this work was developed by Diltz & Böttcher (2014). The calculations are performed in the comoving frame of the emission region, and the calculated spectra and lightcurves are subsequently transformed into the observer’s frame with the Doppler factor δ and the redshift z. At each time-step a power-law distribution of electrons with injection luminosity $L_{\mathrm{inj}}^{\mathrm{e}}$, minimum and maximum Lorentz factor ${\gamma }_{\min }^{\mathrm{e}}$ and ${\gamma }_{\max }^{\mathrm{e}}$, and spectral index s^e is injected into the spherical emission region of radius R, which is pervaded by a tangled magnetic field B. The particle distribution is then evolved self-consistently following a Fokker-Planck equation, considering stochastic acceleration, radiative cooling, and catastrophic losses. The Fokker-Planck equation is solved using a Crank-Nicholson scheme. The stochastic acceleration timescale is parameterized as a multiple η_acc of the escape timescale, which itself is a multiple η_esc of the light-crossing timescale R/c. Any of these timescales is independent of energy implying a “hard-sphere” magnetic field turbulence model for the acceleration term. It should be noted that the acceleration in this case merely acts as a mild re-acceleration of particles. The main acceleration is induced by the injection spectrum, which could originate from a small acceleration region (as in the models of, e.g., Weidinger & Spanier 2015; Chen et al. 2015) that is not accounted for here. The electrons are subject to synchrotron and inverse-Compton cooling including SSC, IC/Disk, IC/BLR, and IC/DT. The latter three depend on the distance r of the emission region from the black hole. The final particle distribution at the end of each time-step is considered in the next time-step with new particles injected on top, and the cycle repeats. It should be noted that the emitted radiation is also self-consistently calculated following the radiative transport equation. This implies that not all emitted photons leave the emission region instantaneously in each time step. Some remain behind and are used in the next time step for all mentioned processes. Eventually an equilibrium solution is found for the particles, where injection, acceleration, and losses balance. Subsequently, any parameter may be disturbed for one or more time steps in order to produce an outburst, after which the code follows the particle evolution until the original equilibrium solution is reached again.

The code has been slightly expanded to include the absorption of γ-rays while they traverse the external photon fields. This adopts the methodology of Böttcher & Els (2016). Additionally, the BLR and DT routines have been slightly expanded to allow for anisotropic photon distributions.

Several constraints can be inferred from the data. The high observed luminosity of the flare along with the short variability implies a large Doppler factor δ in order to keep the particle energy densities low. Unfortunately, no direct constraint on the value of the Doppler factor is available. However, observations of moving radio knots revealed apparent speeds of up to ∼21c (Lister et al. 2013) in the radio jet of 3C 279, also implying large Doppler factors. For the (main) modeling, δ = 30 is adopted, which is well within bounds of the observed apparent superluminal motion (see also Hayashida et al. 2015).

The standard constraint on the size of the emission region is by equating the characteristic flare timescale with the light-crossing timescale of the emission region. Using the value of the characteristic flare timescale from Sect. 3.4, the radius R becomes

$$\begin{array}{c}\hfill R\le \frac{\delta t_{\mathrm{char}}c}{1+z}=1.8\times {10}^{16}\left(\frac{\delta }{30}\right)\mathrm{cm}·\end{array}$$(A.1)

This is the maximum value allowed by the characteristic timescale.

The spectral index of the electron distribution is directly related to the spectral index of the synchrotron component. In the strong cooling regime, the spectral index, α, and the electron spectral index s^e are related by s^e = 2 − 2α, where the electron distribution is n^e(γ) ∝ γ^−se. It has been verified a posteriori that the cooling is indeed in the strong cooling domain. Using the spectral index values for the IR to UV regime from Table 3, the electron spectral index during the Preflare period is s^e = 2.94 ± 0.01, while during the flare it is s^e = 3.12 ± 0.03. For the latter, the average value of Night 1 and Night 2 has been used, since they are consistent within errors. The electron indices are softer than expected by conventional acceleration processes. However, they are in line with typical electron indices derived for 3C 279 (Böttcher et al. 2009) and other FSRQs (e.g., Vercellone et al. 2011; Barnacka et al. 2014; Zacharias et al. 2019).

Simple considerations of the IC process, especially with thermal target photons, lead to the estimate that in a restricted energy range the resulting IC component depends similarly on the electron spectral index as in the synchrotron component. Hence, the spectral index in the γ-ray domain probed by Fermi-LAT, α = 2 − Γ_LAT, should be comparable to the spectral index in the IR to UV domain. Tables 1 and 3 indicate that for the Preflare period the indices are similar, while during the flare the hardening in the γ-ray domain does not correspond to the softening in the IR to UV regime. This can be mitigated by increasing the minimum electron Lorentz factor ${\gamma }_{\min }^{\mathrm{e}}$ during the flare.

In order to model the variability, four parameters have been varied as given in Eqs. (4)–(7). The changes are inspired by either direct measurements, as in the case of the spectral index, or inferences from the changes in the spectrum, such as the Compton dominance. The latter implies a change in the ratio from the external photon density to the magnetic field density. It is assumed that the external thermal fields do not change during the short flare. Hence, the magnetic field must decrease to account for an increase in the Compton dominance. As the Fermi-LAT spectra are almost identical in Night 1 and 2, the Compton dominance is the same in these two nights, which is why the same magnetic field strength is required in both nights (giving the reason for, in total, three Heaviside functions in Eq. (4)). The requirement to increase ${\gamma }_{\min }^{\mathrm{e}}$ has been mentioned before. The flux changes in the synchrotron component imply an increase in particle energy density in order to compensate the decrease in the magnetic field.

Figures A.1 and A.2 show additional leptonic models with a larger Doppler factor, and a larger distance from the black hole, respectively. The respective parameter sets are given in Tables A.1 and A.2. The variability follows the same dependencies as given in Eqs. (4)–(7).

Fig. A.1. Leptonic model using a larger Doppler factor along with the parameters in Table A.1. Data and model lines as in Fig. 7.

Fig. A.2. Leptonic model using a larger distance from the black hole and a larger Doppler factor. The parameters are given in Table A.2. Data and model lines as in Fig. 7.

Appendix B: Lepto-hadronic code and model constraints

The time-dependent lepto-hadronic code used in this work was developed by Diltz & Böttcher (2016) and extended to include external photon fields in Zacharias et al. (2019). This includes the same possibilities as in the leptonic code. Namely, the absorption of γ-rays in the external fields and anisotropic external fields. The code works similarly to the leptonic code described above with the addition of the proton distribution and related effects. In addition to the electrons, protons are injected at each time step with a power-law distribution with injection luminosity $L_{\mathrm{inj}}^{\mathrm{p}}$, minimum and maximum Lorentz factor ${\gamma }_{\min }^{\mathrm{p}}$ and ${\gamma }_{\max }^{\mathrm{p}}$, and spectral index s^p. The protons follow the Fokker-Planck equation with the same structure as the electrons. However, next to synchrotron cooling, protons can also interact with ambient photon fields to produce pions. While the neutral pions are assumed to instantaneously decay into γ-rays, the charged pions decay into muons, which subsequently decay into electrons or positrons. The evolution of the charged pions and muons is also calculated by a Fokker-Planck equation, considering the same effects as for the protons and electrons. The electrons and positrons from the muon decay are used as an additional injection term for the electron evolution. All charged particles are subject to radiative cooling, which is considered to be synchrotron emission. For electrons Compton losses on the ambient fields are also considered. It turns out that these are subdominant. The neutrino spectra produced during pion and muon decay are also calculated. The time-dependency of the code is achieved as in the leptonic case through variations of a few parameters.

Below, the constraints for the lepto-hadronic model are described. Several of the leptonic constraints are reused, most notably the Doppler factor and the size of the emission region.

The spectral indices for the particle distributions can be derived from the observed spectral indices of the observed spectrum listed in Table 3. In fact, for the electrons the result is unchanged. From the interpolated X-ray to γ-ray spectrum one can deduce the proton spectral index assuming slow cooling of the protons. This has also been verified a posteriori. The relation between the observed spectral index α and the proton spectral index s^p is s^p = 3 − 2α, where the proton distribution is described by n^p(γ) = γ^−sp. Using the values of α listed in Table 3, the proton spectral index for Night 1 is 2.16 ± 0.04, 2.10 ± 0.02 for the Maximum, and 2.11 ± 0.04 for Night 2. These are compatible within errors, and are kept constant during the modeling.

The hardening of the HE γ-ray spectrum is mimicked by increasing the maximum proton Lorentz factor during the flare. The apparent break in the X-ray domain is accounted for by a large minimum proton Lorentz factor.

Figure B.1 shows an additional lepto-hadronic model with a larger Doppler factor. The parameter set is given in Table B.1. The variability follows the same dependencies as given in Eqs. (8)–(11).

Fig. B.1. Lepto-hadronic model using a larger Doppler factor along with the parameters in Table B.1. Data and model lines as in Fig. 8.

All Tables

All Figures

	Fig. 1. Observed H.E.S.S. photon spectra for six data sets as labeled. Arrows mark upper limits (99% confidence level). The gray butterfly is the 1σ statistical uncertainty band of the 2015/Night 2 data set. Error bars are statistical only. The second label gives the telescope participation and the analysis used.
In the text

Fig. 2. Observed multiwavelength lightcurves. (a) H.E.S.S. lightcurve derived above an energy threshold of 200 GeV in night-wise time bins with array configuration as indicated. Arrows mark upper limits (99% confidence level). (b) Fermi-LAT lightcurve integrated above 100 MeV in 3 h bins. Gray arrows mark upper limits (95% confidence level). (c) HE γ-ray photon index measured with Fermi-LAT in 3 h bins. (d) Swift-XRT lightcurve integrated between 2 and 10 keV for individual pointings. (e) Optical R band lightcurve from ATOM and SMARTS for individual pointings. (f) Spectral index between the J and B band using SMARTS observations for individual pointings. In all panels, only statistical error bars are shown.

In the text

Fig. 3. Observed multiwavelength lightcurves zoomed in on Night 2. (a) H.E.S.S. lightcurve derived above an energy threshold of 100 GeV in run-wise time bins. (b) Fermi-LAT lightcurve integrated above 100 MeV in 3 h bins. (c) HE γ-ray photon index measured with Fermi-LAT in 3 h bins. (d) Swift-XRT lightcurve integrated between 2 and 10 keV for individual pointings. (e) Optical R band lightcurve from ATOM and SMARTS for individual pointings. In all panels, only statistical error bars are shown, while horizontal bars mark the observation time.

In the text

Fig. 4. Observed multiwavelength SED for the considered time frames with black dots for the Preflare time frame, red filled squares for Night 1, green open squares for the Maximum, and blue diamonds for Night 2. The γ-ray data have been corrected for EBL absorption using the model by Franceschini et al. (2008). The solid lines show a power-law interpolation for the X-ray to γ-ray spectrum, as described in the text.

In the text

Fig. 5. Optical depths for γ-rays emitted along the jet axis at different distances r interacting with photons of the BLR emission lines. Left: BLR modeled with the shell geometry. The crossing of lines at low energies is due to numerical inaccuracies. Right: BLR modeled with the ring geometry. The structure in the optical depth are caused by the contributions of different emission lines to the overall optical depth.

In the text

Fig. 6. Left: best-fit spectra for the BLR ring geometry to the combined Fermi-LAT and H.E.S.S. data. Both data sets are corrected for EBL absorption following Franceschini et al. (2008). The spectral shapes do not change significantly, if the shell geometry is assumed instead. Right: likelihood profile as a function of the distance r for the different assumed intrinsic spectra and BLR geometries.

In the text

Fig. 7. Multiwavelength spectra and models for the four time frames: Preflare (black dots), Night 1 (red filled squares), Maximum (green open squares), and Night 2 (blue diamonds). The γ-ray data have been corrected for EBL absorption using the model by Franceschini et al. (2008). The thick solid lines mark the leptonic models. The thin lines mark spectral components for the Preflare period as labeled.

In the text

Fig. 8. Multiwavelength spectra and models for the four considered time frames: Preflare (black dots), Night 1 (red filled squares), Maximum (green open squares), and Night 2 (blue diamonds). The γ-ray data have been corrected for EBL absorption using the model by Franceschini et al. (2008). The thick solid lines mark the lepto-hadronic photon models, while the thick dashed lines mark the neutrino models. The thin lines mark the spectral components for the Preflare period as labeled.

In the text

	Fig. 9. Lower limits on the sub-luminal linear (left) and quadratic (right) terms of the modified dispersion relation obtained with H.E.S.S. for several AGN as a function of redshift. 1: PKS 2155−304 (H.E.S.S. Collaboration 2011), 2: PG 1553+113 (H.E.S.S. Collaboration 2015), 3: Mrk 501 (H.E.S.S. Collaboration 2019), 4: 3C 279 (this work).
In the text

	Fig. A.1. Leptonic model using a larger Doppler factor along with the parameters in Table A.1. Data and model lines as in Fig. 7.
In the text

	Fig. A.2. Leptonic model using a larger distance from the black hole and a larger Doppler factor. The parameters are given in Table A.2. Data and model lines as in Fig. 7.
In the text

	Fig. B.1. Lepto-hadronic model using a larger Doppler factor along with the parameters in Table B.1. Data and model lines as in Fig. 8.
In the text