© The Authors 2024

1. Introduction

Blazars, a class of jetted active galactic nuclei (AGN) with a relativistic jet aligned very closely with the line of sight to the Earth, are ideal targets for studying the extreme physics of relativistic outflows. These objects emit radiation from the radio band up to the γ-ray range, with highly variable flux (up to a factor ∼10 or even ∼100; e.g., Ackermann et al. 2016) across the entire electromagnetic spectrum on timescales ranging from as short as minutes and hours (e.g., Albert et al. 2007; Hayashida et al. 2015) to as long as months and even years (e.g., Zacharias et al. 2017). The short timescale (t_var ≲ 1 week) flux variability, generally referred to as flaring activity, is of particular interest as it enables us to probe the physical processes occurring in blazar jets at their most extreme.

Blazars are categorised as BL Lacertae objects (BL Lacs) and flat-spectrum radio quasars (FSRQs) based on the prominence of emission lines in their optical spectra. FSRQs, with noticeable optical emission lines, generally have higher luminosities primarily dominated by γ-ray emission. The spectral energy distributions (SEDs) of blazars typically display two bumps: a low-energy synchrotron radiation bump, and a high-energy bump, the nature of which remains unclear. In the leptonic scenario, this bump is explained as being due to the inverse Compton (IC) scattering of low-energy photons by the same electron population that produces the synchrotron emission, with the low-energy photons being either synchrotron photons (i.e., the synchrotron self-Compton (SSC) scenario; e.g., Maraschi et al. 1992) or external photons from an accretion disk (e.g., Dermer et al. 1992), a broad line region (BLR; e.g., Sikora et al. 1994), and/or a dusty torus (e.g., Błażejowski et al. 2000) (i.e., the external Compton (EC) scenario). The SED of BL Lac objects is better explained by the SSC scenario, while the emission from FSRQs is more in line with the EC scenario (e.g., Tavecchio et al. 1998; Abdo et al. 2010a; Böttcher et al. 2013). Additionally, hadronic scenarios have been proposed, in which the γ-ray emission of blazars is generated via processes involving hadrons (e.g., Mücke & Protheroe 2001; Mücke et al. 2003).

The origin of the blazar flaring behaviour remains a mystery despite numerous studies, and several mechanisms have been proposed to explain this phenomenon. One of the simplest descriptions is a leptonic one-zone model, in which the broad-band emission of blazars originates from a compact region in the jet (a ‘blob’) filled with electron–positron plasma and moving relativistically along the jet axis (e.g., Katarzyński et al. 2001). Flares arise from perturbations of this configuration, either due to microphysics inside the emitting zone or the macrophysical properties of the blob, such as its geometry or kinematics. In the first case, a flare can be produced due to various transient processes such as particle injection (e.g., Albert et al. 2007) and/or acceleration due to shocks (e.g., Marscher & Gear 1985; Sikora et al. 2001; Böttcher & Baring 2019), turbulence (e.g., Tammi & Duffy 2009; Asano et al. 2014; Baring et al. 2017), or magnetic reconnection (e.g., Giannios et al. 2009; Petropoulou et al. 2016; Shukla & Mannheim 2020). In the second case, emission can be enhanced due to, for example, an increase in the particle number density or in the Doppler factor of the blob (e.g., Casadio et al. 2015; Paliya et al. 2015; Larionov et al. 2016; Luashvili et al. 2023), with the latter caused by either the increase in the bulk Lorentz factor of the blob (e.g., Ghisellini & Tavecchio 2008) or changes in the viewing angle when the emitting zone moves along a twisted or bent jet (Abdo et al. 2010b; Raiteri et al. 2017) or a helical trajectory (Villata & Raiteri 1999).

Understanding blazar flaring behaviour involves modelling the broad-band emission during flaring states. Two alternative approaches exist to describe the evolution of the emitting particle population (leptonic and/or hadronic) due to various physical processes and/or changes of the physical conditions in the source: the kinetic (e.g., Kardashev 1962; Kirk et al. 1998) or Monte-Carlo (MC) approach (e.g., Mücke et al. 1999; Summerlin & Baring 2012), or the combination of the two (e.g., Chen et al. 2011; Dimitrakoudis et al. 2012). Various numerical codes developed for blazar modelling implement these approaches. In particular, numerical codes based on the kinetic approach track the evolution of the spectrum of the particle population residing in the emitting region (and in some models also in the accelerating region) (e.g., Mastichiadis & Kirk 1995; Chiaberge & Ghisellini 1999; Katarzyński et al. 2003; Dmytriiev et al. 2021) or in the large-scale blazar jet (e.g., Zacharias et al. 2022), and calculate the time-dependent particle spectrum evolving due to various physical processes, as well as the associated varying multi-wavelength (MWL) emission. In the present work, we specifically focus on the kinetic description of blazar variability within the leptonic framework.

Particles producing the synchrotron and γ-ray emission lose energy, which is referred to as particle cooling. Within the leptonic scenario, electrons in the blazar emitting zone suffer energy losses via two mechanisms: (1) synchrotron cooling, and (2) IC processes, in which electrons transfer part of their energy to a photon, upscattering it to high energies (IC cooling). The second mechanism becomes especially important during extreme γ-ray flares. When the IC process proceeds in the Thomson regime (γx ≪ 1, where γ is the electron Lorentz factor and is the dimensionless energy of the seed photon), electrons lose only a small fraction of their energy, ΔE/E ≪ 1, and the cooling can be well approximated as a continuous process in the kinetic description. However, when the IC scattering proceeds in the Klein–Nishina (KN) regime (γx ∼ 1), electrons lose a substantial fraction of their energy in a single interaction, and the process cannot be treated as continuous.

The full transport equation for the latter case has an integro-differential form (e.g., Blumenthal & Gould 1970), which presents an increased challenge compared to equations in partial derivatives. To remedy this problem, several authors have derived continuous-loss approximations in an attempt to include KN effects, providing a reasonable description of the IC cooling (e.g., Böttcher et al. 1997; Moderski et al. 2005). As a result, the non-continuous effects of cooling have previously been neglected in blazar emission models. This was commonly accepted because Zdziarski (1988) showed that substantial differences in the electron spectrum shape only appear in the case of mono-energetic injection and target photon field spectra. Consequently, most existing blazar emission modelling codes employ a continuous-loss term in the kinetic equation to describe IC cooling, despite its validity being limited to cases with small relative energy losses per scattering event. This is a noteworthy simplification, given that blazars often exhibit significant Klein–Nishina effects.

In this study, we revisit the full integro-differential equation for particle IC cooling in the KN regime, which properly treats the large relative jumps in energy experienced by the electrons, and investigate deviations from the continuous-loss approximation in the context of blazars. We concentrate on scenarios where discrete cooling effects are most pronounced. Specifically, we examine the flaring states of FSRQs, which are characterised by strong Compton dominance (U_rad ≫ U_B), and consider the target photon field of the broad line region (BLR), whose spectral shape closely resembles a monoenergetic distribution. We identify distinctive features in the electron distribution and SED induced by the discrete jumps in energy experienced by the electrons, and assess the overall significance of such effects in blazars.

2. Transport equation and numerical approach

In this section, we examine the non-continuous cooling terms in the kinetic equation, and present the numerical implementation for solving it.

2.1. Transport equation with the full cooling term

In our model, we assume a one-zone leptonic emission scenario, with electron-positron plasma in the compact region of the jet (a ‘blob’) producing the observed varying broad-band emission. We refer to both electrons and positrons as electrons hereafter. The γ-ray emission in this scenario is produced via the IC process and the electrons therefore experience IC cooling. A standard form of the kinetic (Fokker-Planck) equation describing the evolution of the electron spectrum N_e(γ, t) in the blazar emitting zone due to particle injection, escape, and cooling – with the latter treated as a continuous process –, is given by (e.g., Chiaberge & Ghisellini 1999)

(1)

where t_esc is the characteristic timescale of particle escape, Q_inj(γ, t) is the spectrum of injected particles, and the quantity is the continuous cooling rate, which is a sum of synchrotron and IC cooling rates, and , respectively:

(2)

(3)

where U_B = B²/(8π) is the magnetic field density, σ_T is the Thomson cross-section, m_e is the electron rest mass, and c is the velocity of light in vacuum. For the IC cooling rate , a useful continuous-loss approximation was derived by Moderski et al. (2005):

(4)

with

(5)

where represents the energy density of the target photon field (per unit of photon energy interval), with both energy density and target photon energy x being in the frame of the blob, and x_min/max being minimum and maximum dimensionless energies of soft photons, respectively. The above-mentioned approximation attempts to reasonably treat the KN effects, in particular, it correctly takes into account the decrease in the cross-section with increasing seed photon energy. However, it still fails to capture large jumps in energy experienced by the particles due to the inherent non-continuous nature of this effect. The full transport equation to treat large jumps in energy experienced by particles in an exact manner may be written as (e.g., Blumenthal & Gould 1970; Zdziarski 1988)

(6)

with

(7)

being the Compton kernel by Jones (1968), which expresses the rate (per unit of time) of electron downscattering from γ to γ′, with being the number density of the target photons (per unit of photon energy interval) in the frame of the blob, and

(8)

(9)

(10)

The first integral term on the right-hand side of Eq. (6) represents IC downscattering from γ to lower Lorentz factors, while the second integral term describes IC downscattering from higher Lorentz factors to γ. The synchrotron losses are still represented by a continuous-loss term , which is reasonable given the (always) small fractional energy loss in this process.

2.2. Numerical implementation

To solve the full transport equation (Eq. (6)), we use and extend the existing numerical code EMBLEM by Dmytriiev et al. (2021), which solves the standard kinetic equation in Eq. (1) (with additional Fermi-I and Fermi-II acceleration terms) using the Chang & Cooper (1970) numerical scheme, which is well-suited to the above-mentioned equation form in particular, with energy losses given by a continuous term. A number of modifications were made to the EMBLEM code in order to adapt it for solving the equation with a full cooling term.

The Compton kernel (Eq. (7)) is sharply peaked around the region γ ≈ γ′, which makes numerical integration challenging. We follow the approach proposed by Zdziarski (1988), where the integration domain spanning from γ = 1 to γ → ∞ is split into three subregions, with the middle peculiar region, γ/(1 + δ)≤γ′≤γ(1 + δ), δ ≪ 1, being treated separately. By performing a Taylor expansion around the peculiar point γ = γ′, the term for the middle region can be reduced to a continuous-loss form (of the same form as the first term on the right-hand side of Eq. (1)). Following this approach, the full IC cooling term is then expressed via three different terms (Zdziarski 1988):

(11)

with representing the IC cooling transition rate in the continuous-loss regime:

(12)

It is worth noting that the rate is not the same as (Eq. (4)), as the latter was designed as an attempt to approximate the combined effect provided by all three terms of Eq. (11) in a continuous framework, while the rate only incorporates small energy losses suffered by electrons. To avoid troublesome numerical integration with a sharply peaked integrand, the Eq. (12) is integrated analytically over the Lorentz factors (Zdziarski 1988), which yields

(13)

where s = 4xγ and g = min(δ/s, 1).

To adapt the EMBLEM code to the integro-differential transport equation (Eq. (6)), we first transform it into a form where the terms have the same mathematical structure as the standard kinetic equation (Eq. (1)) that the original code is designed to solve, relying on the Chang & Cooper (1970) numerical scheme. This is achieved by numerically computing the non-continuous cooling terms in Eq. (11) on the Lorentz factor grid, at each time step, substituting the sought quantity N_e(γ, t) with the electron spectrum from the previous time step (ensuring a sufficiently small time step on the time grid). The result represents a certain tabulated function of γ and can therefore be treated as a time-dependent source (injection-like) term in addition to Q_inj(γ). For the continuous part of the full cooling term (Eq. (11)), we compute the IC cooling rate in the continuous regime using Eq. (13). The integral over the seed photon energies is evaluated numerically and the result is added to the synchrotron cooling rate (Eq. (2)) in the code, yielding the total continuous cooling rate, treated in the same manner as in Eq. (1). As a result, we simplify the integro-differential equation to a standard differential one, enabling the application of the Chang & Cooper (1970) scheme. After solving the equation, we obtain the updated electron spectrum, which is then used to reevaluate the non-continuous cooling terms on the Lorentz factor grid, improving the accuracy of the estimated values. This process is repeated iteratively until convergence is achieved. Throughout our numerical computations, we use δ = 0.05. For each solution, we ensure that the normalizations (total number of particles) are equal between the electron spectra for non-continuous and continuous-loss cooling scenarios.

3. Application to the FSRQ 3C 279

In this section, we apply the code developed above to simulate the electron spectrum and SED during a realistic flaring state, representing a recent flare of the FSRQ 3C 279, and compare the results from a continuous cooling case with those from a non-continuous cooling case.

3.1. 3C 279: The archetypal FSRQ

3C 279 (z = 0.5362) is a prototypical FSRQ that exhibits intense and complex variability across all wavebands, including the γ-ray range. It is one of the most studied blazars thanks to numerous MWL campaigns on the target (Hayashida et al. 2012, 2015). The multi-frequency SED of 3C 279 displays a characteristic two-bump structure with the low-energy component peaking in the infrared band, and the high-energy component reaching its maximum in the domain ranging from ∼0.1 GeV to a few GeV, depending on the spectral state. 3C 279 has been continuously observed by the Fermi Gamma-Ray Space Telescope since 2008. During a high state in 2013–2014, the source showed a very hard γ-ray spectrum with a high level of Compton dominance of ∼300, and a peak γ-ray flux of f_γ(> 0.1 GeV)≃10⁻⁵ ph cm⁻² s⁻¹. In June 2015, it showed an extreme flare with the historically highest peak in γ-ray flux, f_γ(> 0.1 GeV)≃(3.6 ± 0.2)×10⁻⁵ ph cm⁻² s⁻¹, and with a very short flux-doubling timescale of ≲5 min (Ackermann et al. 2016).

3.2. Physical scenario for the low and flaring states

3.2.1. General model

For our study, we chose the most extreme flare of 3C 279, specifically that observed in June 2015 with the historically highest γ-ray flux of F_γ(0.1 − 100 GeV)∼2 × 10⁻⁸ erg cm⁻² s⁻¹, as well as the highest ever Compton dominance level of ∼800 (Fig. 1 in Dmytriiev et al. 2023). We adhere to the one-zone leptonic scenario, assuming that the varying emission of the source is produced within a spherical ‘blob’ of fixed radius R_b that is situated at a distance r_emz from the central black hole and is travelling down the jet at a velocity close to the speed of light, characterised by a bulk Lorentz factor Γ and a Doppler factor δ_b. To reduce the number of free parameters, we assume the viewing angle of the jet i_j such that δ_b = Γ at all times. The γ-rays are produced in the blob via IC scattering of a target photon field. This target photon field comprises two components: an external field of the BLR (time-independent) and the synchrotron radiation field (which is time dependent). In our case of FSRQs, the external photon field dominates the high-energy emission production and IC cooling processes. We model the BLR photon field simply as a single Gaussian-shaped Ly α emission line with a fixed central energy of E_Lyα, c = 10.2 eV and fixed width ΔE_Lyα = 0.2 eV corresponding to an average velocity of BLR clouds of ∼6000 km s⁻¹ (both energies given in the BLR frame). The luminosity of the BLR is considered to be a constant fraction ξ_BLR of the accretion disk luminosity, for which we adopt the value L_AD = 6 × 10⁴⁵ erg s⁻¹ (Hayashida et al. 2015). The calculation of the energy density of the external photon field in the frame of the blob follows the approach by Hayashida et al. (2012). We assume that particles are injected into the blob with a log-parabola electron injection spectrum above a certain minimum injection Lorentz factor γ_inj, c,

(14)

where γ₀ = 200 is the (arbitrarily chosen) pivot Lorentz factor. We attribute the lower-energy cutoff γ_inj, c in the injection spectrum to the electron–proton co-acceleration occurring in the shocked region of the jet (refer to e.g., Zech & Lemoine 2021); however, we do not model this process in the present study. It is worth mentioning that we also attempted the modelling using a power-law injection spectrum with and without an exponential cutoff. However, we were not able to reproduce the low-state SED data with such injection spectra. Finally, particles also escape the emitting region at a characteristic timescale of t_esc ∼ 1 R_b/c. Throughout this work, we use a redshift of 3C 279 of z = 0.5362, and the Hubble parameter value of H₀ = 70 km s⁻¹ Mpc⁻¹.

3.2.2. Low-state model

Following this framework, we first modelled the low state of the source using the multi-band data set from Hayashida et al. (2012). We simulate the low state as an asymptotic steady-state established due to the competition of injection, cooling, and particle escape (based on the approach by e.g., Böttcher et al. 2013; Dmytriiev et al. 2021). In this case, softening of the electron spectrum due to cooling arises naturally without the need to artificially introduce this effect. In addition, in our modelling, we disregard the low-energy radio measurements, assuming that emission from the source in this regime is dominated by the extended jet rather than the compact emitting zone.

3.2.3. June 2015 flare model

Next, we performed a modelling of the June 2015 flaring state. We used the MWL data collected during the dedicated observational campaign, which include spectral measurements in the optical-UV (UVOT instrument), X-ray (Swift-XRT), γ-ray (Fermi-LAT), and very high energy (VHE) γ-ray (H.E.S.S.) bands reported by H. E. S. S. Collaboration (2019). Our primary focus is to accurately reproduce the SED at the flare peak (the ‘Maximum’ time frame in H. E. S. S. Collaboration 2019, MJD 57189.125 – 57189.25), while also attempting to provide a satisfactory fit of the SED ∼8.9 h prior to the peak at the early stage of the outburst (the ‘Night 1’ time frame in H. E. S. S. Collaboration 2019, MJD 57188.756 – 57188.88) within the self-consistent time-dependent framework. However, our modelling approach does not include fitting the shapes of MWL light curves. Despite this, our methodology ensures a reasonable approximation of the dynamic timescales around the flare peak.

Following the framework introduced by Dmytriiev et al. (2021), we treat the flaring state as a perturbation of the low-state configuration induced by a transient physical process (e.g., particle acceleration, injection, etc.), which may be accompanied by a change in one or a few global physical parameters. To minimise complexity, we adhere to the one-zone scenario with transient Fermi-II reacceleration of the particle population within the blob due to intervening turbulence simulated using a dedicated setup of the EMBLEM code (see Sect. 6.1 of Dmytriiev et al. 2021). We vary the Fermi-II acceleration timescale t_acc, FII and the escape timescale t_esc, FII, with the latter typically longer than the low-state escape timescale due to particle diffusion through the turbulent magnetic field (e.g., Tramacere et al. 2011). Another parameter is the total duration of the perturbative acceleration episode t_dur, FII, which is constrained based on the time span of the γ-ray flux rise observed in the Fermi-LAT light curve of the flare (H. E. S. S. Collaboration 2019), of namely ∼12 h.

We also assessed the significance of the triplet pair production (TPP) process (e⁻ + γ → e⁻ + e⁻ + e⁺) under our conditions. At high electron and/or target photon energies, TPP losses start competing with IC cooling losses (e.g., Mastichiadis et al. 1994). Using developments by Mastichiadis et al. (1994), we determine that for electron interactions with BLR target photons, boosted in the emitting zone frame with a Doppler factor of δ_b ∼ 10–100, TPP energy losses () dominate over IC losses at γ > 10⁹. As no electrons with Lorentz factors exceeding ∼10⁷ exist in the emitting region, we conclude that the TPP loss process has a negligible effect on the electron spectrum in our case.

3.3. Results

Initially, we used the continuous-loss approximation for cooling by Moderski et al. (2005) to reproduce the observed characteristics of both the quiescent and flaring state. Our model successfully replicates the low-state data set, as illustrated in Fig. 1, with the corresponding model parameters listed in Table 1. The underlying electron spectrum is depicted in Fig. 1 (blue dashed curves).

Fig. 1. Comparison of the low-state SEDs and electron spectra for the two cooling descriptions. The blue dashed curves represent the model with continuous-loss approximation for cooling (Moderski et al. 2005), whereas green solid curves indicate the same model, but with the full non-continuous cooling term. Left panel: modelling of the low-state SED of 3C 279. The black points display the low-state data set of the source (Hayashida et al. 2012). The adjacent narrow bottom panel provides a zoom onto the SED ratio of the non-continuous to continuous-loss cooling description (in linear scale). Right panel: underlying low-state electron spectra for the two cooling descriptions. The electron spectra are shown in γ3Ne(γ) representation to better highlight the discrepancies between the two cooling models. Adjacent narrow bottom panels show the ratio of the electron spectrum between the non-continuous and continuous-loss cooling scenarios (in linear scale), with the second lower panel showing a zoom onto the Lorentz factor domain above 100.

For the flaring state, initially, we attempted to replicate the MWL SED data of the flare peak by adjusting t_acc, FII and t_esc, FII, but found that the model consistently underpredicts the observed γ-ray flux level and Compton dominance. To address this, we also varied the Doppler factor δ_b, fl (with the bulk Lorentz factor again adhering to Γ_fl = δ_b, fl) and magnetic field B_fl during the acceleration episode, given that Γ_fl and B_fl control the Compton dominance, , with U_rad being the density of the BLR field in the AGN rest frame. The bulk Lorentz factor may be enhanced due to passage of the blob through an active part of the jet (e.g., Ghisellini & Tavecchio 2008), with the Doppler factor increasing to reach δ_b, fl = Γ_fl owing to a modest decrease in the jet viewing angle because of, for example, the curvature of this jet segment. The magnetic field decrease may arise from magnetic reconnection and/or turbulent dissipation of magnetic energy in this jet region. Furthermore, a study by Dmytriiev et al. (2023) concludes that the states of extreme Compton dominance in 3C 279 are most likely explained by strong variations in the magnetic field and moderate variations in the bulk Lorentz factor. The (total) perturbation duration in the blob frame is well constrained and is adjusted only in a narrow range t_dur, FII = (12 ± 1) h × δ_b, fl/(1 + z).

As a result, we find t_acc, FII = 3.6 R_b/c, t_esc, FII = 4.2 R_b/c, δ_b, fl = 40, and B_fl = 0.2 G are necessary to adequately fit the flare peak data, as well as the pre-peak state (see Table 2). The acceleration and escape timescales appear to be in overall agreement with the doubling and halving timescales observed in the Fermi-LAT light curve. The ‘Night 1’ (pre-peak) dataset is fit by the resulting (time-dependent) model at the time moment t₁ = 3.6 R_b/c (time elapsed after the acceleration onset), and the dataset during the peak with the same model at t₂ = 12.8 R_b/c, with the time difference t₂ − t₁ correctly matching the ∼8.9 h (observer’s frame) interval between the two observations. The best-fit (total) duration of the perturbation, t_dur, FII = t₂ = 12.8 R_b/c, translates to ≈12.3 h in the observer’s frame, and is well compatible with the observed flux rise time. We therefore connected the low state, the state in the initial stage of the flare, and the state during the flare peak in a single model within a self-consistent time-dependent framework. The modelling results for the two high states are depicted in Fig. 2. Our model offers a generally reasonable description of the observed data. However, while it reproduces the overall γ-ray emission level, there is a notable underprediction of the flux below 1 GeV. This discrepancy likely arises from our simplified treatment of the BLR photon field, which only considers a single emission line instead of the actual broader spectrum. Additionally, a contribution from the dusty torus field may be necessary, particularly for γ-ray production at lower energies, as suggested by previous studies (e.g., Hayashida et al. 2012; Dermer et al. 2014).

Fig. 2. Modelling of the June 2015 flaring state of 3C 279. The black, blue, and purple data points display the MWL data during the low state, flare peak, and in the pre-peak state of ‘Night 1’ (≈8.9 h earlier), respectively (H. E. S. S. Collaboration 2019). The data points at the highest energies are uncorrected for the extragalactic background light (EBL) absorption. The error bars for the optical and X-ray data points are relatively small and are not visible. The grey curve shows the low-state model from Fig. 1, the cyan curve indicates the time-dependent flare model at the moment t1 = 3.6 Rb/c, fitting the pre-peak state data, and the green curve illustrates the same model at the moment t2 = 12.8 Rb/c, fitting the flare peak data. The model SEDs were absorbed on the EBL using the model by Domínguez et al. (2011). The model curves are displayed in the dashed style to highlight the use of the continuous-loss approximation in the model calculation.

It is worth noting that, while blazar emission models often feature parameter degeneracies, our specific model exhibits relatively independent parameters. In particular, the Doppler factor regulates both synchrotron and γ-ray peak fluxes, with a different scaling relation (which may be somewhat complex due to strong cooling effects), while the magnetic field directly controls the synchrotron peak only. The acceleration timescale influences the resulting spectral index (‘hardness’) in a given spectral band, and the escape timescale affects the total number of particles in the system, N_tot. This number has a direct impact on the (relative) level of X-ray flux, which in our scenario is produced mostly by the SSC process, for which the emission is , rather than ∝N_tot as for synchrotron or external Compton processes. This approach yields a seemingly unique solution within our chosen physical scenario.

Comparing with alternative models for this flare, H. E. S. S. Collaboration (2019) model the peak SED within the leptonic external Compton model by a step-change in four parameters, namely the magnetic field, electron injection luminosity, minimum electron Lorentz factor, and the index of the electron spectrum. Importantly, our scenario, while also incorporating a decrease in the magnetic field, allows the electron spectrum to naturally evolve due to the Fermi-II reacceleration process, providing a more self-consistent approach.

Subsequent to this, we repeated all the above-mentioned simulations of both the low and flaring state using the derived parameters, this time incorporating the full (non-continuous) IC cooling description. A comparison of the underlying electron spectra and SEDs between the non-continuous and continuous-loss scenarios for cooling for the low state is presented in Fig. 1, and the same, this time for the flaring state, is shown in Fig. 3.

Fig. 3. Comparison of the electron spectrum (left) and the associated SED (right) between two cooling descriptions for the peak of the June 2015 flaring state of 3C 279. The electron spectra here are again shown in γ3Ne(γ) representation. The bottom panel of the electron spectrum–SED plot displays the ratio between the electron spectra–SEDs in the case of non-continuous cooling (orange solid curves) and the continuous-loss approximation (green dashed curves) in linear scale.

4. Discussion

Comparing the two cooling scenarios in the low state allows us to reveal steady-state implications of non-continuous effects in a configuration that results from a balance of physical processes. In contrast, the flaring state comparison uncovers dynamic effects arising from transient changes in source conditions and interplay between various physical processes.

4.1. Low state

The comparison for the low state (Fig. 1) reveals general agreement in electron spectra for Lorentz factors γ ≳ 100. Discrepancies are below ∼10%, with a moderate pile-up observed around the Lorentz factor corresponding to the KN transition, γ ≈ 2000 (Fig. 1). In the non-continuous cooling scenario (solid curves), the electron spectrum appears harder in the injection-driven regime, showing a ‘turning point’ at γ ≈ 600, aligning closely with the anticipated position of the cooling break. Substantial effects emerge below γ_inj, c = 130, where no injection takes place, including a prominent ‘tail’ of particles below γ ∼ 10 and a particle excess above the tail (10 ≲ γ ≲ 130).

Collectively, the discrepancies in the electron spectra result in differences in the observed SED of up to ≲10% (top panels of Fig. 1). The most prominent features in the SED for the non-continuous cooling case (compared to continuous-loss case) include a hardening of the spectral slope in the γ-ray band beyond the peak of the high-energy bump.

The hardening of the electron spectrum around the cooling break, , together with the excess immediately above around the KN transition, γ_KNx ∼ 1, can be attributed to the fact that the continuous-loss approximation might provide an inadequate description of particle cooling in the KN transition regime, tending to overestimate the cooling effect. This inaccuracy, in turn, might arise from the discrete nature of particle cooling. At the same time, this approximation overlooks large relative jumps from the region around γ_KN to significantly lower Lorentz factors; a fraction of particles experience these jumps while interacting in the KN regime. This results in an excess and a ‘tail’ feature at very low energies. We conclude therefore that this tail feature might represent an ‘echo’ of the particle cooling close to the KN transition regime.

4.2. Flaring state

An immediate observation is the substantial increase in discrepancies between the two cooling scenarios in the flaring state compared to the low state. We note that the low-energy tail and narrow excess observed in the low-state electron spectrum are now replaced by an extended ‘shelf’-like excess by a (roughly) constant factor of ∼2, stretching over the Lorentz factors below γ ∼ 100 (Fig. 3). Below the expected Lorentz factor of the KN transition (γ_KN ≈ 1300), a spectral softening is evident, transitioning to spectral hardening above this point. Subsequently, an extreme and relatively narrow pile-up, reaching a factor of ∼10, emerges around γ ∼ 2 × 10⁴. Finally, at the highest Lorentz factors, we observe deviations by a factor of ∼1.3–2.

These deviations induce strong differences in the SED across the entire frequency domain, with broad-range deviations of up to ≃50%, along with two narrow pile-ups by factors of ∼3 and ∼8 (Fig. 3). The most substantial discrepancies are observed in the UV, X-ray, and VHE bands, along with additional notable disparities in the GeV γ-ray band. Similarly to the low state, the SED in the non-continuous cooling case exhibits a harder slope beyond the peak of the high-energy bump (GeV-to-VHE regime). Remarkably, the non-continuous cooling model predicts, on average, a several-fold higher flux in the VHE regime during the flare. Upon examination of the SEDs, it is evident that the model SED for the non-continuous cooling case can no longer adequately describe the data. Specifically, the model SED for the non-continuous cooling case exhibits a markedly softer spectral slope in the X-rays, with discrepancies surpassing the (relatively small) uncertainties in the available X-ray data by a very large margin. Important deviations are also observed in the optical–UV band. In the GeV γ-ray band, discrepancies exceed the Fermi-LAT uncertainties, particularly in the lower energy range of approximately 1–5 GeV. The non-continuous cooling effects therefore appear to be significant.

Now we proceed to an interpretation of the observed features and discrepancies. The substantially enhanced discrepancies between the two cooling scenarios primarily arise from an increase in the overall number of particles in the flaring state (as compared to the low state), as well as an increase in the relative number of particles at higher Lorentz factors, particularly beyond the Klein–Nishina threshold (γ > γ_KN), as a result of the reacceleration process. Additionally, a slightly increased Lorentz and Doppler factor during the flaring state leads to a boost of the seed photon field in the blob frame and amplification of the cooling rate, as well a boost in the target photon energy. Consequently, the majority of the γ-ray production now occurs much closer to the KN regime, enhancing the discrete cooling effects.

For the domain below γ_inj, c, unlike in the low state, it is now no longer dominated by losses only (cooling and escape). Instead, Fermi-II acceleration now competes with cooling effects. With the inferred Fermi-II timescale, 3.6 R_b/c, acceleration dominates for γ ≤ 50, that is, t_{cool, syn + IC}(γ)> t_acc, FII, with . Hence, the shelf-like feature below γ ∼ 100 can be attributed to the Fermi-II process reaccelerating low-energy particles downscattered from interactions in the KN regime, effectively spreading them into an extended homogeneous excess.

Above the KN transition, as the IC cooling loses its efficiency, the total cooling timescale increases with a higher γ and eventually becomes comparable to the acceleration timescale again (t_{cool, syn + IC} ∼ t_acc, FII), occurring around γ ≈ 3 × 10⁴, where the cooling timescale exhibits a local maximum. Beyond this point, synchrotron cooling begins to dominate the total cooling rate, with the cooling timescale decreasing again with an increasing γ, resulting in synchrotron cooling strongly prevailing above the acceleration process (γ > 10⁵). In the Lorentz factor range where cooling is the most inefficient (γ ∼ 10⁴–10⁵), a distinct ‘bump’ appears in the electron spectra due to the maximised relative importance of reacceleration. The position of this bump for the non-continuous cooling scenario is appreciably shifted towards higher Lorentz factors. This can be explained by the continuous-loss approximation potentially underestimating the decline in the efficiency of the IC process in the KN regime, which leads to overestimation of the overall cooling effect. The shift in the position of the bump leads to the observed strong pile-up at γ ≈ 3 × 10⁴, as well as the pronounced narrow pile-up appearing in the low- and high-energy SED component. At the highest Lorentz factors, the electron spectrum for the non-continuous cooling case consistently appears higher (showing an extended particle excess), which is likely also due to the overestimation of the cooling effect in the continuous-loss scenario.

5. Conclusions

In this work, we carried out a novel and detailed exploration of the non-continuous (discrete) IC cooling effects in γ-ray(Compton)-dominated blazars (FSRQs) arising due to IC interactions proceeding in the KN regime. Our analysis centres on the archetypal FSRQ 3C 279, and involves modelling of the brightest γ-ray flare of this source observed in June 2015. We find that the discrete cooling effects may become important and significantly affect the electron spectra and SEDs of blazars, inducing a range of distinguishable features, such as low-energy tails, narrow and extended excesses and pile-ups, spectral softening and hardening, shifts in the position of the cooling bump or break, and so on. More specifically, our study has unveiled the following main effects:

1. Low state: In the injection-dominated regime, the electron spectra agree well, with discrepancies below 10%. The most notable features emerge at low Lorentz factors in the injection-free regime, and include a significant particle excess and a low-energy ‘tail’, which appear to be low-energy reflections of the KN transition. We interpret these features as due to particles interacting in the KN regime and experiencing large jumps to much lower Lorentz factors, an effect that is not properly taken into account by the continuous-loss description. These deviations lead to minor SED differences of up to 10%. Therefore, non-continuous cooling effects can be neglected in the low states of blazars.

2. Flaring state: During flares, the electron spectrum and SED show significantly amplified discrepancies between continuous and non-continuous cooling, yielding a distinct array of features. Specifically, we investigated a (typical) flare scenario, where the flux increase is initiated by transient particle reacceleration (turbulent or Fermi-II mechanism in this case), with acceleration competing with cooling. In the flaring state, a pronounced extended particle excess emerges at low Lorentz factors, stemming from particles downscattered to very low Lorentz factors and subsequently reaccelerated. Furthermore, a moderate excess at high Lorentz factors, as well as a significant shift in the position of the cooling bump or break (resulting from a competition between acceleration and different regimes of cooling), suggest an overestimation of cooling effects in the KN regime by the continuous-loss approximation. Broad-range SED deviations peak at 50%, with the most prominent ones seen in the hard X-ray and VHE bands. Also, the shape of the high-energy SED bump (γ-ray domain) displays a noticeable distortion, including a considerable hardening of the SED in the GeV-to-VHE domain in the non-continuous case, as well as a strong, narrow pile-up. The discrepancies between the two cooling scenarios are found to exceed the uncertainties of the observational data by a considerable margin, notably in the X-ray domain. This highlights the importance of the non-continuous effects during flares and demonstrates that their imprint on the emission spectra can be effectively discerned in the modelling of data from current MWL instruments; we speculate that these imprints will be even more easily discernable in the data provided by future-generation observatories like CTA.

In summary, taking into account the non-continuous cooling effects is essential for accurate modelling of blazar flaring states, as these effects become non-negligible during strong FSRQ flares characterised by high Compton dominance.

Acknowledgments

The work of M.B. was supported by the South African Research Chairs Initiative of the National Research Foundation (any opinion, finding, and conclusion or recommendation expressed in this material is that of the authors, and the NRF does not accept any liability in this regard) and the Department of Science and Innovation of South Africa through SARChI grant no. 64789. A.D. thanks Andrzej Zdziarski for useful discussions on the subject. A.D. acknowledges the FSK-GAMMA-1 computer cluster facility at the Centre for Space Research (North-West University) which was particularly helpful for performing computationally-expensive simulations for this project. Also, A.D. thanks Pieter Van der Merwe and Robert Brose for their ideas and suggestions that were of great help for the numerical implementation of the full transport equation and for the code optimization. Finally, the authors extend their gratitude to the anonymous referee for their critical and thorough review of the manuscript, as well as for providing important requests and suggestions that significantly enhanced the quality of the presented work.

References

All Tables

All Figures

Fig. 1. Comparison of the low-state SEDs and electron spectra for the two cooling descriptions. The blue dashed curves represent the model with continuous-loss approximation for cooling (Moderski et al. 2005), whereas green solid curves indicate the same model, but with the full non-continuous cooling term. Left panel: modelling of the low-state SED of 3C 279. The black points display the low-state data set of the source (Hayashida et al. 2012). The adjacent narrow bottom panel provides a zoom onto the SED ratio of the non-continuous to continuous-loss cooling description (in linear scale). Right panel: underlying low-state electron spectra for the two cooling descriptions. The electron spectra are shown in γ3Ne(γ) representation to better highlight the discrepancies between the two cooling models. Adjacent narrow bottom panels show the ratio of the electron spectrum between the non-continuous and continuous-loss cooling scenarios (in linear scale), with the second lower panel showing a zoom onto the Lorentz factor domain above 100.

In the text

Fig. 2. Modelling of the June 2015 flaring state of 3C 279. The black, blue, and purple data points display the MWL data during the low state, flare peak, and in the pre-peak state of ‘Night 1’ (≈8.9 h earlier), respectively (H. E. S. S. Collaboration 2019). The data points at the highest energies are uncorrected for the extragalactic background light (EBL) absorption. The error bars for the optical and X-ray data points are relatively small and are not visible. The grey curve shows the low-state model from Fig. 1, the cyan curve indicates the time-dependent flare model at the moment t1 = 3.6 Rb/c, fitting the pre-peak state data, and the green curve illustrates the same model at the moment t2 = 12.8 Rb/c, fitting the flare peak data. The model SEDs were absorbed on the EBL using the model by Domínguez et al. (2011). The model curves are displayed in the dashed style to highlight the use of the continuous-loss approximation in the model calculation.

In the text

Fig. 3. Comparison of the electron spectrum (left) and the associated SED (right) between two cooling descriptions for the peak of the June 2015 flaring state of 3C 279. The electron spectra here are again shown in γ3Ne(γ) representation. The bottom panel of the electron spectrum–SED plot displays the ratio between the electron spectra–SEDs in the case of non-continuous cooling (orange solid curves) and the continuous-loss approximation (green dashed curves) in linear scale.

In the text