A Steady-state Spectral Model for Electron Acceleration and Cooling in Blazar Jets: Application to 3C 279

The broadband $\nu {F}_{\nu }$ spectral energy distributions (SEDs) of blazars have two broad components ("bumps"): a low-frequency one peaking between the IR and X-rays, and a high-frequency one peaking in the $\gamma$-rays. The lower-frequency component is almost certainly synchrotron emission from nonthermal electrons (including positrons) accelerated in the jet. The radio emission in blazars is thought to be generated outside of the primary (γ-ray) emission region, due to the superposition of synchrotron spectra from a range of jet locations where the optical depth ${\tau }_{\mathrm{opt}}\approx 1$ (e.g., Blandford & Königl 1979; Konigl 1981), and cannot be formed by the same component as the low-frequency peak emission process because the former must be self-absorbed. The origin of the higher energy component is a matter of debate. It could be due to Compton scattering of synchrotron photons by the same nonthermal electrons that produce the synchrotron (synchrotron self-Compton, or SSC), which is thought to produce the $\gamma$-rays in high-peaked BL Lac objects (e.g., Maraschi et al. 1992; Dermer & Schlickeiser 1993; Sikora et al. 1994; Bloom & Marscher 1996; Błażejowski et al. 2000). Alternatively, it could be due to the Compton scattering of lower energy seed photons impinging on the jet (external Compton or EC) from an accretion disk (e.g., Dermer et al. 1992; Dermer & Schlickeiser 1993), broad-line region (BLR; e.g., Sikora et al. 1994; Blandford & Levinson 1995; Ghisellini & Madau 1996), or dust torus (e.g., Kataoka et al. 1999; Błażejowski et al. 2000). The seed photon source is closely tied to the location in the jet of the $\gamma$-ray-emitting region. For an emitting region at distances $r\lesssim 0.01\ \mathrm{pc}$ from the black hole (BH), the scattering of disk photons dominates; for an emitting region in the range $0.01\lesssim r\lesssim 0.1\ \mathrm{pc}$, the scattering of BLR photons dominates; and for distances $0.1\lesssim r\lesssim 1\ \mathrm{pc}$, the scattering of dust torus photons dominates.

Alternatively, the $\gamma$-rays could be produced by hadronic processes, such as proton synchrotron (e.g., Aharonian 2000; Mücke & Protheroe 2001; Mücke et al. 2003; Reimer et al. 2004) or by the decay products of proton–photon interactions (e.g., Sikora et al. 1987; Mannheim & Biermann 1992; Mannheim 1993; Protheroe 1995) with protons co-accelerated in the jet with the electrons, and the synchrotron emission of the resulting particles produced by these hadronic processes. Hadronic processes for $\gamma$-ray production are in some cases disfavored due to excessive energetics (e.g., Böttcher et al. 2013; Zdziarski et al. 2015; Petropoulou & Dermer 2016), leaving electron Compton scattering as the most likely source for $\gamma$-rays.

The electron acceleration mechanism operating during the observed blazar flares is currently not well understood (e.g., Madejski & Sikora 2016; Romero et al. 2017). The emitting electrons are probably energized by some combination of shock acceleration (e.g., Summerlin & Baring 2012; Marscher 2014), stochastic scattering (e.g., Katarzyński et al. 2006; Lefa et al. 2011; Asano & Hayashida 2015; Baring et al. 2017), or electrostatic acceleration due to magnetic reconnection (e.g., Giannios et al. 2009; Giannios 2013; Petropoulou et al. 2016; Sironi et al. 2016). Often, rather than modeling the acceleration mechanism, a functional form for the emitting electron distribution is assumed, and theoretical SEDs are calculated from this (e.g., Tavecchio et al. 1998; Bednarek & Protheroe 1997, 1999; Finke et al. 2008; Dermer et al. 2009; Hayashida et al. 2012).

We endeavor to improve on these types of models by determining the electron distribution using a self-consistent framework, in which we treat particle acceleration, escape, and energy losses by solving a Fokker–Planck equation. We include diffusive shock and stochastic acceleration, as well as adiabatic and radiative losses. In particular, including the full Compton cross-section is more physically realistic and can be important in the case where the high-energy SED peak is brighter compared with the low-energy peak. The Klein–Nishina (KN) turnover can cause inefficient losses at high energies, leading to a harder electron distribution and hence harder resulting synchrotron and Compton spectra (e.g., Dermer & Atoyan 2002; Moderski et al. 2005). Dermer et al. (2014) studied the use of a log-parabola electron distribution for modeling blazars, in particular for modeling the SEDs of the FSRQ 3C 279 in an equipartition framework. These authors had difficulty explaining the $\gamma$-ray emission detected from the SEDs by the Fermi Large Area Telescope (LAT) at $\gtrsim 5\ \mathrm{GeV}$. All of the models we present here are able to explain the $\geqslant 5\,\mathrm{GeV}$ $\gamma$-rays detected by the Fermi-LAT. Other authors (e.g., Asano et al. 2014; Diltz & Böttcher 2014; Asano & Hayashida 2015) have performed time-dependent modeling of several blazar flares, including second-order particle acceleration and energy losses described by the full Compton cross-section. We have created a similar model, although we use a time-independent, steady-state analysis for simplicity. However, we include a more complete physical modeling of the BLR (Finke 2016) and scattering of dust torus photons, both being important seed photon sources for $\gamma$-ray production.

Dermer et al. (2014) applied their modeling technique to the well-studied FSRQ 3C 279. This was the first blazar detected by EGRET on board the Compton Gamma-Ray Observatory in 1991 (Hartman et al. 1992) and only the second blazar detected in $\gamma$-rays after 3C 273 by COS B (Swanenburg et al. 1978). In 2008, Fermi-LAT began collecting data for 3C 279 as part of its standard survey mode, providing regular monitoring of this source and reinvigorating multiwavelength efforts. The first multiwavelength campaign for 3C 279 found a large optical polarization swing associated with a $\gamma$-ray outburst detected by Fermi-LAT (Abdo et al. 2010). Hayashida et al. (2012) performed further analysis on several spectral epochs. They modeled the broadband SEDs with broken power-law electron distributions and two scenarios: a propagating emission region model and a model with two distinct synchrotron emission regions to account for the high- and low-energy features simultaneously. The sequence of SEDs provided by Hayashida et al. (2012) is relatively complete, in both wavelength and time, and provides an excellent resource for SED modeling.

The FSRQ 3C 279 has a redshift z = 0.536, giving it a luminosity distance ${d}_{L}=3100\ \mathrm{Mpc}$ in a cosmology where ${H}_{0}=70\ \mathrm{km}\,{{\rm{s}}}^{-1}\ {\mathrm{Mpc}}^{-1}$, ${{\rm{\Omega }}}_{m}=0.3$, and ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.7$. We use this distance in all of our relevant calculations.

Our goal in this paper is to develop a new one-zone leptonic model, in which the jet emission is assumed to come predominantly from one emitting region—a homogeneous spherical "blob." We therefore neglect the radio emission, which is produced elsewhere (e.g., Blandford & Königl 1979; Konigl 1981). The propagation angle relative to the line of sight, θ, introduces a Doppler factor ${\delta }_{{\rm{D}}}={[{\rm{\Gamma }}(1-\beta \cos \theta )]}^{-1}$. We will assume here that ${\delta }_{{\rm{D}}}={\rm{\Gamma }}$. The blob contains nonthermal electrons, randomly oriented (isotropic) in the frame co-moving with the blob (the "co-moving frame"). We assume that the co-moving blob radius ${R}_{b}^{{\prime} }$ is constrained by the light-crossing timescale t_v, so that ${R}_{b}^{{\prime} }\lesssim c{\delta }_{{\rm{D}}}{t}_{v}/(1+z)$. The blob moves at a highly relativistic speed, $v=\beta c$, giving it a bulk Lorentz factor ${\rm{\Gamma }}={(1-{\beta }^{2})}^{-1/2}$. The electrons are accelerated by diffusive shock and stochastic processes inside the blob, which also contains a tangled, turbulent magnetic field. The blob may also contain protons, but these are radiatively unimportant in our model.

The blob is embedded in a jet launched from a supermassive BH at the center of the active galaxy. Matter accretes onto the BH, creating a thermal accretion disk with temperature profile and emission spectrum given by Shakura & Sunyaev (1973). The disk's emission sometimes appears as a "blue bump" in the optical/UV spectra of quasars. We will refer to the frame that is stationary with respect to the BH as the "stationary frame."

The optical spectra of active galactic nuclei often contain broad emission lines, indicating the presence of highly ionized, fast-moving gas in a BLR. Reverberation mapping indicates the BLR gas is traveling in Keplerian orbits around the BH, and different lines emit primarily at different distances from the BH (e.g., Peterson & Wandel 1999). The BLR absorbs disk photons and re-emits them as line radiation (the eponymous broad lines). Some of these BLR photons make their way to the jet blob, where they are Compton-scattered by the nonthermal electrons there. The exact geometry of the BLR is not well known. In blazar spectral modeling, it is sometimes assumed that Lyα will be the dominant emission source from the BLR, and that other emission lines can be ignored (e.g., Moderski et al. 2003; Dermer et al. 2014) or that the whole BLR spectrum can be represented by a thermal distribution at a single radius (e.g., Ghisellini et al. 1998; Ghisellini & Tavecchio 2008; Tavecchio & Ghisellini 2008).

Joshi et al. (2014) used 35 lines in their BLR model, assuming all of the line emission was produced throughout a single shell regardless of line. In contrast, Finke (2016) created a stratified BLR model, with 26 different lines emitting at different radii, in a spherical shell or flattened configuration with the BH at the center (see also Abolmasov & Poutanen 2017), and we adopt the same picture here. The radii of the different line sources were based on an SDSS composite quasar spectrum (Vanden Berk et al. 2001). Finke (2016) provided formulae for computing the accurate Compton-scattering spectra from these lines, taking into account the complicated geometry of the BLR. We model the stratification of the BLR in a simplified fashion, using an isotropic approximation for the Compton-scattering formulae, and allow the energy density from a particular monochromatic line to decrease outside of the radius where that line is emitted. Thus, in this picture, the most prominent emission line for Compton scattering, out of the 26 considered here, depends on the location of the blob.

We also include a dust torus in our model, which absorbs and re-emits a fraction of the disk photons. The infrared emission from the torus is modeled similarly to a BLR line as a radiation field isotropic in the stationary frame that is monochromatic, with photon energy given by the peak of a blackbody spectrum. Some subset of the infrared dust photons become seed photons for Compton scattering by electrons in the blob.

Our model is based on a quantitative description of the basic physical processes experienced by the electrons contained in a magnetically confined plasma blob in the emitting region of the blazar jet. We are especially interested in formulating a self-consistent correspondence between the underlying electron distribution and the observed spectral distribution, where practicable. Specifically, the electrons experience first- and second-order Fermi acceleration, Bohm diffusion (energy-dependent escape), adiabatic losses, synchrotron losses, and EC losses due to interactions with seed radiation from the disk, dust torus, and the BLR. The spectral model also includes synchrotron self-absorption (SSA) and SSC processes. SSC losses were neglected in the solutions to the electron transport equation since this process renders the equation nonlinear (e.g., Zacharias & Schlickeiser 2012a, 2012b; Zacharias 2014). However, these losses are unlikely to have a strong effect, since EC losses will dominate. We consider the full Compton cross-section in our transport equation solution (e.g., Böttcher et al. 1997) and spectral calculations (e.g., Jones 1968; Blumenthal & Gould 1970), where the KN regime becomes important at higher energies.

We compare our steady-state model to multiwavelength SED data for 3C 279, originally published in Hayashida et al. (2012). Some of the light curves for these SEDs were first published in Abdo et al. (2010), but not all of the SEDs we use were included, and they did not model the SEDs in that paper. Hayashida et al. (2012) use a broken power-law or double broken power-law electron distribution as necessary, with both one- and two-zone models. They used the X-ray data collected during flares as upper limits (because the X-ray and γ-ray light curves during the flare were not correlated) when they attempted to compare their model with the multiwavelength SED. Four of the SEDs from Hayashida et al. (2012) were also analyzed by Dermer et al. (2014). The latter used a log-parabola electron distribution to minimize the number of free model parameters. They also assumed equipartition between the energy densities of the magnetic field, photons, and electrons, which minimizes the jet power, and approximated the BLR contribution to EC as Lyα, all in an attempt to ground the SED model interpretation in physics. After all of that, the $\gamma$-ray data disagree with the Dermer et al. (2014) model above  5 GeV, prompting the authors to suggest a second synchrotron or hadronic component as possible remedies.

We build on each of these models in analyzing the same spectral data (2008 August 4–2009 February 23), but we calculate our electron distribution by solving the electron transport equation in a more physically realistic, fully self-consistent model. We do not assume equipartition between any energy densities and have a more physically realistic, stratified BLR model.

The details of our model calculation can be found in Section 2, including the transport equation electron distribution solution. The spectral calculation methods for each individual component are in Section 3. In Section 4, we perform several parameter studies to illustrate the capabilities of the electron distribution solution for physical interpretation and highlight the flexibility of the spectral calculations. We present our comparison to 3C 279 in Section 5 and discuss the significance of our results in Section 6.

The evolution of electrons in the jet blob can be described by a Fokker–Planck equation. We use this equation to treat particle energization via diffusive shock and stochastic acceleration, energy losses from adiabatic and radiative processes, particle escape via Bohm diffusion, and the continuous injection of monoenergetic electrons. In this subsection, all of our calculations occur in the co-moving frame, which we generally denote by primes (e.g., ${N}_{e}^{{\prime} }({\gamma }^{{\prime} })$) but suppress for simplicity when there are no other frames under consideration. Our Fokker–Planck equation takes the form

$$\begin{eqnarray}\displaystyle \frac{\partial {N}_{e}}{\partial t} & = & \displaystyle \frac{{\partial }^{2}}{\partial {\gamma }^{2}}\left(\displaystyle \frac{1}{2}\displaystyle \frac{d{\sigma }^{2}}{{dt}}{N}_{e}\right)-\displaystyle \frac{\partial }{\partial \gamma }\left(\left\langle \displaystyle \frac{d\gamma }{{dt}}\right\rangle {N}_{e}\right)\\ & & -\displaystyle \frac{{N}_{e}\gamma {D}_{0}}{\tau }+{\dot{N}}_{\mathrm{inj}}\delta (\gamma -{\gamma }_{\mathrm{inj}}),\end{eqnarray} \tag{ 1 }$$

where ${N}_{e}(\gamma ,t)$ is the electron number distribution and $\gamma \equiv E/({m}_{e}{c}^{2})$ is the electron Lorentz factor, with m_e and c denoting the electron mass and the speed of light, respectively. The quantities ${\dot{N}}_{\mathrm{inj}}$ and ${\gamma }_{\mathrm{inj}}$ in Equation (1) represent the particle injection rate and the Lorentz factor of the injected monoenergetic electrons, respectively.

We can obtain the solution to the steady-state problem by setting

$$\begin{eqnarray}&&\displaystyle \frac{\partial {N}_{e}}{\partial t}=0,\end{eqnarray} \tag{ 2 }$$

so that ${N}_{e}(\gamma ,t)\to {N}_{e}(\gamma )$. We describe the individual terms below, while detailed derivations can be found in Lewis et al. (2016).

The particle injection mechanism is not well understood, but we assume here that the electrons accelerated inside the blob originate as members of a low-energy thermal distribution and that the acceleration and emission regions are co-spatial. Since the model considered here includes a significant component of stochastic particle acceleration, represented by the second-order term in Equation (1), the steady-state electron distribution (the solution to Equation (2)) will lose any "memory" of the energy of the injected electrons. In this scenario, the steady-state energy distribution is independent of the precise form of the energy distribution of the injected electrons. We can take advantage of this insensitivity by assuming that the injected electrons have a monoenergetic distribution with a very low injection energy. The monoenergetic injection is represented by the δ-function in Equation (1).

We adopt the hard-sphere approximation, so that the Fokker–Planck "broadening coefficient" in Equation (1) is given by

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}\displaystyle \frac{d{\sigma }^{2}}{{dt}}={D}_{0}\,{\gamma }^{2},\end{eqnarray} \tag{ 3 }$$

where ${D}_{0}\propto {{\rm{s}}}^{-1}$ is the energy-diffusion coefficient (Park & Petrosian 1995). We use two forms of the drift coefficient, describing the mean net acceleration rate. They differ in their description of the Compton-cooling term. In the Thomson approximation, the drift coefficient takes the form

$$\begin{eqnarray}&&\left\langle \displaystyle \frac{d\gamma }{{dt}}\right\rangle ={D}_{0}\left[4\gamma +a\gamma -{b}_{\mathrm{syn}}{\gamma }^{2}-{\gamma }^{2}\displaystyle \sum _{j=1}^{J}{b}_{{\rm{C}}}^{(j)}\right],\end{eqnarray} \tag{ 4 }$$

where D₀ is the hard-sphere scattering coefficient for collisions with MHD waves. The first term on the right-hand side of Equation (4) expresses the mean energization rate due to stochastic acceleration, and the term containing a accounts for first-order Fermi acceleration at the shock as well as adiabatic cooling. Shock acceleration dominates over adiabatic losses if $a\gt 0;$ otherwise, adiabatic cooling dominates.

Synchrotron cooling is expressed by the term containing ${b}_{\mathrm{syn}}$ in Equation (4), which is defined by

$$\begin{eqnarray}&&{b}_{\mathrm{syn}}\equiv \displaystyle \frac{| \,{\dot{{\gamma }}}_{\mathrm{syn}}| }{{D}_{0}{\gamma }^{2}}=\displaystyle \frac{{\sigma }_{{\rm{T}}}{B}^{2}}{6\pi {m}_{e}{{cD}}_{0}},\end{eqnarray} \tag{ 5 }$$

where ${\sigma }_{{\rm{T}}}$ is the Thomson cross-section. The term containing ${b}_{{\rm{C}}}^{(j)}$ in Equation (4) expresses energy losses due to the scattering of external seed photons impinging on the jet and colliding with the electrons in the blob. In our calculation, there are multiple sources of seed photons for Compton scattering, and each is considered independently. Hence, the Compton term is a sum over J different seed photon sources, with associated dimensionless coefficients

$$\begin{eqnarray}&&{b}_{{\rm{C}}}^{(j)}\equiv \displaystyle \frac{| {\dot{\gamma }}_{{\rm{C}}}| }{{D}_{0}{\gamma }^{2}}=\displaystyle \frac{4{\sigma }_{{\rm{T}}}{{\rm{\Gamma }}}^{2}{u}_{\mathrm{ph}}^{(j)}}{3{m}_{e}{{cD}}_{0}},\end{eqnarray} \tag{ 6 }$$

where ${u}_{\mathrm{ph}}^{(j)}$ is the photon energy density for each seed photon source (see Section 3.3.6).

We also consider a more advanced version of the model in which we utilize the full Compton cross-section, including the KN correction. In this case, the expression for the drift coefficient can be written as

$$\begin{eqnarray}&&\left\langle \displaystyle \frac{d\gamma }{{dt}}\right\rangle ={D}_{0}\left[4\gamma +a\gamma -{b}_{\mathrm{syn}}{\gamma }^{2}-{\gamma }^{2}\displaystyle \sum _{j=1}^{J}{b}_{{\rm{C}}}^{(j)}H(\gamma {\epsilon }_{\mathrm{ph}}^{(j)})\right],\end{eqnarray} \tag{ 7 }$$

where ${\epsilon }_{\mathrm{ph}}^{(j)}$ is the incident photon energy in units of m_ec². The function H gives the energy-loss rate with the full Compton cross-section included and is defined by

$$\begin{eqnarray}&&H(y)\equiv \displaystyle \frac{9}{32}\displaystyle \frac{1}{{y}^{2}}{G}_{\mathrm{BMS}}(y),\end{eqnarray} \tag{ 8 }$$

where (Böttcher et al. 1997)

$$\begin{eqnarray}{G}_{\mathrm{BMS}}(y) & = & \displaystyle \frac{8}{3}y\displaystyle \frac{1+5y}{{(1+4y)}^{2}}-\displaystyle \frac{4y}{1+4y}\left(\displaystyle \frac{2}{3}+\displaystyle \frac{1}{2y}+\displaystyle \frac{1}{8{y}^{2}}\right)\\ & & +\mathrm{ln}(1+4y)\left[1+\displaystyle \frac{3}{y}+\displaystyle \frac{3}{4{y}^{2}}+\displaystyle \frac{\mathrm{ln}(1+4y)}{2y}\right.\\ & & \left.-\displaystyle \frac{\mathrm{ln}(4y)}{y}\right]\\ & & -\displaystyle \frac{5}{2y}-\displaystyle \frac{{\pi }^{2}}{6y}-2+\displaystyle \frac{1}{y}\displaystyle \sum _{n\,=\,1}^{\infty }\displaystyle \frac{{(1+4y)}^{-n}}{{n}^{2}}.\end{eqnarray} \tag{ 9 }$$

We note that $H(y)\to 1$ in the Thomson limit $y\ll 1$ (see Figure 1). Thus, Equation (7) reduces to Equation (4) in the Thomson regime.

Zoom In Zoom Out Reset image size

Following Lewis et al. (2016), our transport equation also includes a term describing the energy-dependent escape of electrons from the blob via Bohm diffusion, which can be written as

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial {N}_{e}}{\partial t}\right|}_{\mathrm{esc}}=-\displaystyle \frac{{N}_{e}\gamma {D}_{0}}{\tau },\end{eqnarray} \tag{ 10 }$$

where the dimensionless escape timescale constant τ is defined in the Bohm limit as

$$\begin{eqnarray}&&\tau \equiv \displaystyle \frac{{R}_{b}^{{\prime} 2}{{qBD}}_{0}}{{m}_{e}{c}^{3}}\ \end{eqnarray} \tag{ 11 }$$

(Lewis et al. 2016). The final term on the right-hand side of Equation (1) is the source term, representing the injection of ${\dot{N}}_{\mathrm{inj}}$ electrons per second with Lorentz factor ${\gamma }_{\mathrm{inj}}$. The value of ${\dot{N}}_{\mathrm{inj}}$ is related to the electron injection luminosity, ${L}_{{\rm{e}},\mathrm{inj}}$, by

$$\begin{eqnarray}&&{L}_{{\rm{e}},\mathrm{inj}}={m}_{e}{c}^{2}{\gamma }_{\mathrm{inj}}{\dot{N}}_{\mathrm{inj}}.\end{eqnarray} \tag{ 12 }$$

The Compton-cooling term in Equation (1) is an addition to the drift coefficient terms previously treated in Lewis et al. (2016). While the previous equation was analytically solvable, the inclusion of the KN formalism leaves the transport equation without an analytical solution, motivating the numerical treatment employed here (see Section 2.3).

In Lewis et al. (2016), we have previously solved Equation (1) analytically in the Thomson regime ($\gamma {\epsilon }_{\mathrm{ph}}\ll 1$), i.e., if the drift coefficient is given by Equation (4). Our solution, stated in the context of this paper, can be written as

$$\begin{eqnarray}{N}_{e}(\gamma ) & = & \displaystyle \frac{{\dot{N}}_{\mathrm{inj}}{e}^{{b}_{\mathrm{tot}}{\gamma }_{\mathrm{inj}}/2}}{{b}_{\mathrm{tot}}{D}_{0}{\gamma }_{\mathrm{inj}}^{2+a/2}}\displaystyle \frac{{\rm{\Gamma }}(\sigma -\lambda +1/2)}{{\rm{\Gamma }}(1+2\sigma )}{e}^{-{b}_{\mathrm{tot}}\gamma /2}{\gamma }^{a/2}\\ & & \times \,\left\{\begin{array}{l}{M}_{\lambda ,\sigma }({b}_{\mathrm{tot}}\gamma ){W}_{\lambda ,\sigma }({b}_{\mathrm{tot}}{\gamma }_{\mathrm{inj}}),\gamma \leqslant {\gamma }_{\mathrm{inj}}\\ {M}_{\lambda ,\sigma }({b}_{\mathrm{tot}}{\gamma }_{\mathrm{inj}}){W}_{\lambda ,\sigma }({b}_{\mathrm{tot}}\gamma ),\gamma \geqslant {\gamma }_{\mathrm{inj}}\end{array}\right.\end{eqnarray} \tag{ 13 }$$

(Lewis et al. 2016), where the total loss coefficient, ${b}_{\mathrm{tot}}$, is defined by

$$\begin{eqnarray}&&{b}_{\mathrm{tot}}\equiv {b}_{\mathrm{syn}}+\displaystyle \sum _{j=1}^{J}{b}_{{\rm{C}}}^{(j)}.\end{eqnarray} \tag{ 14 }$$

The indices of the Whittaker functions, ${M}_{\lambda ,\sigma }$ and ${W}_{\lambda ,\sigma }$, in the branching solution are given by (e.g., Abramowitz & Stegun 1972)

$$\begin{eqnarray}&&\lambda =2-\displaystyle \frac{1}{{b}_{\mathrm{tot}}\tau }+\displaystyle \frac{a}{2},\quad \mathrm{and}\ \quad \sigma =\displaystyle \frac{a\,+\,3}{2}.\end{eqnarray} \tag{ 15 }$$

While the form of the solution in Equation (13) is not valid once we consider KN effects, it does show that the electron distribution is dependent upon many of the same physical parameters as the emission spectra, especially the magnetic field and the photon energy densities.

With the full Compton energy-loss term included (Equation (7)), we can no longer solve Equation (1) analytically. Hence, we utilized a fourth-order Runge–Kutta finite-differencing scheme, outlined in Press et al. (1992), modified to treat a second-order ODE.

We employ a grid that is uniform in terms of the logarithm of the electron Lorentz factor γ. The minimum Lorentz factor is ${\gamma }_{\min }={10}^{0}$ in all of our calculations, and the Lorentz factor of the injected electrons is set at ${\gamma }_{\mathrm{inj}}=1.01$. One of our goals here is to attempt to explain the flares' SEDs using a first-principles approach, without any requirement for electron pre-acceleration. Hence, we choose a small value for the electron injection energy in order to model injection from the thermal electron distribution, rather than invoking any form of electron pre-acceleration, such as acceleration at shocks, which is assumed in many previous studies. The maximum Lorentz factor ${\gamma }_{\max }$ is set to accommodate the full dynamic range of the solution and is typically between ${10}^{4.5}$ and ${10}^{7.5}$.

Solutions to this class of transport equation follow two distinct branches for $\gamma \leqslant {\gamma }_{\mathrm{inj}}$ and $\gamma \geqslant {\gamma }_{\mathrm{inj}}$. The solution for the electron distribution, ${N}_{e}(\gamma )$, must be continuous at ${\gamma }_{\mathrm{inj}}$, but the derivative $\partial {N}_{e}/\partial \gamma$ is discontinuous, as discussed in detail by Lewis et al. (2016).

We have developed a process for computing the electron distribution in a piecewise fashion. At the high- and low-energy boundaries, we have used the flux-free condition, i.e., the electrons should not cross the boundary at $\gamma ={\gamma }_{\min }=1$ or $\gamma ={\gamma }_{\max }$. Once the fundamental solutions are obtained in the regions above and below the injection energy, one normally utilizes the derivative-jump condition for ${N}_{e}(\gamma )$ at the injection energy ${\gamma }_{\mathrm{inj}}$. However, we found that the accuracy required by this Wronskian method is prohibitive, due to the stiffness of the differential equation, and therefore we employed an alternative method to normalize the global solution for ${N}_{e}(\gamma )$ that is based on the requirement of a steady-state balance between electron injection and escape. This is discussed in detail in the Appendix. It requires much less precision than the Wronskian method and may be useful in other steady-state applications involving stiff ordinary differential equations. We used the analytic Thomson regime solution (Equation (13)) to verify the machinery for the full Compton solution in the Thomson limit.

In Equation (8), recall that $H(y)=1$ in the Thomson limit ($y\ll 1$), but the realization of this limit depends on both the electron Lorentz factor, $\gamma$, and the incident photon energy, ${\epsilon }_{\mathrm{ph}}$, since $y=\gamma {\epsilon }_{\mathrm{ph}}$ in the full Compton calculation. As the incident photon energy ${\epsilon }_{\mathrm{ph}}$ decreases, the Compton-cooling term eventually approaches the Thomson limit (Figure 1).

It is interesting to contrast the results obtained for the electron distribution ${N}_{e}(\gamma )$ under the influence of pure synchrotron cooling versus synchrotron plus Thomson cooling (Section 2.2), or synchrotron plus the full Compton-cooling rate (Section 2.3) due to a single external photon source with various photon energy values ${\epsilon }_{\mathrm{ph}}$. This comparison is carried out in Figure 2. For our choice of parameters (Tables 1 and 2), for low ${\epsilon }_{\mathrm{ph}}$, where the Compton solution approaches the Thomson regime, and in the Thomson regime, Compton cooling dominates (${{\rm{\Gamma }}}^{2}{u}_{\mathrm{ph}}/[{B}^{2}/(8\pi )]\gg 1$). However, as ${\epsilon }_{\mathrm{ph}}$ increases, the Compton cooling at high $\gamma$ becomes less efficient (Figure 1), and therefore synchrotron losses begin to dominate. For ${\epsilon }_{\mathrm{ph}}\leqslant {10}^{-4}$, the solution is indistinguishable from the synchrotron-only case.

Zoom In Zoom Out Reset image size

Table 1.  Physical and Assumed Parameters for 3C 279

Variable (Units)	Symbol	Value
Energy of Lyα	${\epsilon }_{\mathrm{Ly}\alpha }$	$2.0\times {10}^{-5}$
Luminosity Distance (cm)	${d}_{{\rm{L}}}$	$9.61\times {10}^{27}$
Redshift	z	0.536
Dust Efficiency	ξ	0.1
Particle Injection Energy	${\gamma }_{\mathrm{inj}}$	1.01

Download table as:  ASCIITypeset image

Table 2.  Base Parameters for Parameter Study

Variable (Units)	Symbol	Base Value
Variability Timescale (s)	${t}_{\mathrm{var}}$	$2.0\times {10}^{4}$
Magnetic Field (G)	B	1.0
Doppler Factor	${\delta }_{{\rm{D}}}$	25.0
Energy Density of Lyα (erg cm−3)	${u}_{\mathrm{Ly}\alpha }$	$2.5\times {10}^{-4}$
Dust Temperature (K)	${T}_{\mathrm{dust}}$	1400
Disk Luminosity (erg s−1)	${L}_{\mathrm{disk}}$	$1.0\times {10}^{46}$
Second-order Fermi Coefficient	D0	$2.0\times {10}^{-6}$
First-order Fermi Coefficient	a	−3.8
Energy Injection Rate	${L}_{\mathrm{inj}}$	$8.3\times {10}^{29}$

Download table as:  ASCIITypeset image

The disk spectrum is approximated using the Shakura & Sunyaev (1973) formalism, which yields

$$\begin{eqnarray}&&\epsilon {L}^{\mathrm{disk}}(\epsilon )=1.12\ {\left(\displaystyle \frac{\epsilon }{{\epsilon }_{\max }}\right)}^{4/3}{e}^{-\epsilon /{\epsilon }_{\max }}\end{eqnarray} \tag{ 16 }$$

(e.g., Dermer et al. 2014), where $\epsilon$ is the stationary-frame dimensionless energy and ${m}_{e}{c}^{2}{\epsilon }_{\max }=10\,\mathrm{eV}$. The observed $\nu {F}_{\nu }$ spectrum from the disk will then be

$$\begin{eqnarray}&&{f}_{{\epsilon }_{\mathrm{obs}}}^{\mathrm{disk}}=\displaystyle \frac{(1+z){\epsilon }_{\mathrm{obs}}{L}^{\mathrm{disk}}((1+z){\epsilon }_{\mathrm{obs}})}{4\pi {d}_{L}^{2}},\end{eqnarray} \tag{ 17 }$$

where the factors of $(1+z)$ take into account the cosmological redshift.

The spectrum radiated by the dust torus is assumed to be an infrared blackbody, so that

$$\begin{eqnarray}&&\epsilon {L}^{\mathrm{dust}}(\epsilon )=\displaystyle \frac{15{L}^{\mathrm{dust}}}{{\pi }^{4}}\displaystyle \frac{{(\epsilon /{\rm{\Theta }})}^{4}}{\exp (\epsilon /{\rm{\Theta }})-1},\end{eqnarray} \tag{ 18 }$$

where the dimensionless dust temperature ${\rm{\Theta }}={k}_{{\rm{B}}}{T}_{\mathrm{dust}}/({m}_{e}{c}^{2})$ and ${k}_{{\rm{B}}}$ denotes the Boltzmann constant. Hence, the observed dust torus spectrum is given by

$$\begin{eqnarray}&&{f}_{{\epsilon }_{\mathrm{obs}}}^{\mathrm{dust}}=\displaystyle \frac{(1+z){\epsilon }_{\mathrm{obs}}{L}^{\mathrm{dust}}((1+z){\epsilon }_{\mathrm{obs}})}{4\pi {d}_{L}^{2}}\ \end{eqnarray} \tag{ 19 }$$

(e.g., Dermer et al. 2014).

In addition to the disk and torus components, the observed blazar SED also includes emission from the nonthermal electrons in the plasma blob, which is propagating inside the relativistic jet. These components are due to synchrotron, SSC, and EC-scattered emission. Once we have solved the steady-state transport equation (Equation (1)) to determine the electron distribution in the co-moving frame of the blob, ${N}_{e}^{{\prime} }({\gamma }^{{\prime} })$, we can use that solution to calculate the jet emission components due to synchrotron, SSC, EC/Dust, and EC/BLR. Our goal is to use our steady-state model to interpret the SED of 3C 279 observed at different epochs. Each epoch will have a different set of physical parameters and a different electron distribution.

In general, quantities calculated in the co-moving (blob) frame are primed, and those in the stationary (BH) frame or observer (Earth) frame are unprimed. Quantities in the observer frame will be affected by the relativistic Doppler effect (e.g., Rybicki & Lightman 1979) as well as the cosmological redshift. In general, the dimensionless energy in the observer frame is related to that in the co-moving frame via

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{obs}}=\displaystyle \frac{{\delta }_{{\rm{D}}}}{1\,+\,z}\ {\epsilon }^{{\prime} }.\end{eqnarray} \tag{ 20 }$$

3.3.1. Synchrotron Emission

In the observer frame, the $\nu {F}_{\nu }$ synchrotron emission spectrum is given by

$$\begin{eqnarray}&&{f}_{{\epsilon }_{\mathrm{obs}}}^{\mathrm{syn}}=\displaystyle \frac{\sqrt{3}\epsilon ^{\prime} {\delta }_{{\rm{D}}}^{4}{e}^{3}B}{4\pi {{hd}}_{{\rm{L}}}^{2}}{\int }_{1}^{\infty }d{\gamma }^{{\prime} }{N}_{e}^{{\prime} }({\gamma }^{{\prime} })\ R(x),\end{eqnarray} \tag{ 21 }$$

where

$$\begin{eqnarray}&&x=\displaystyle \frac{4\pi \epsilon ^{\prime} {m}_{e}^{2}{c}^{3}}{3{eBh}{\gamma }^{{\prime} 2}},\end{eqnarray} \tag{ 22 }$$

and

$$\begin{eqnarray}&&R(x)\equiv \displaystyle \frac{x}{2}{\int }_{0}^{\pi }d\theta \ \sin \theta \ {\int }_{x/\sin \theta }^{\infty }{dw}\ {K}_{5/3}(w)\end{eqnarray} \tag{ 23 }$$

(Crusius & Schlickeiser 1986). Here, ${K}_{5/3}(w)$ denotes the modified Bessel function of order 5/3. Numerically, we use the approximate expression for R(x) from Finke et al. (2008); other useful approximations are given by Zirakashvili & Aharonian (2007) and Joshi & Böttcher (2011).

We also consider the self-absorption of the synchrotron radiation, which depletes the low-energy portion of the synchrotron spectrum. The self-absorbed spectrum in the observer frame, denoted by ${f}_{{\epsilon }_{\mathrm{obs}}}^{\mathrm{syn},\mathrm{abs}}$, is computed in the slab approximation using

$$\begin{eqnarray}&&{f}_{{\epsilon }_{\mathrm{obs}}}^{\mathrm{syn},\mathrm{abs}}=\displaystyle \frac{1-{e}^{-{\tau }^{\mathrm{SSA}}}}{{\tau }^{\mathrm{SSA}}}{f}_{{\epsilon }_{\mathrm{obs}}}^{\mathrm{syn}}\end{eqnarray} \tag{ 24 }$$

(e.g., Dermer & Menon 2009), where ${f}_{{\epsilon }_{\mathrm{obs}}}^{\mathrm{syn}}$ is computed using Equation (21) and where the SSA optical depth, ${\tau }^{\mathrm{SSA}}$, is given by

$$\begin{eqnarray}&&{\tau }_{\mathrm{SSA}}=\displaystyle \frac{-{R}_{b}^{{\prime} }{c}^{2}}{8\pi \nu {m}_{e}{c}^{2}}{\int }_{1}^{\infty }{P}_{\nu }(\gamma ){\gamma }^{2}\,\displaystyle \frac{d}{d\gamma }\left[\displaystyle \frac{n(\gamma )}{{\gamma }^{2}}\right]d\gamma \end{eqnarray} \tag{ 25 }$$

(e.g., Böttcher et al. 1997), where

$$\begin{eqnarray}&&{P}_{\nu }(\gamma )=\displaystyle \frac{\sqrt{3}{e}^{3}B}{{m}_{e}{c}^{2}}\,R(x)\end{eqnarray} \tag{ 26 }$$

(e.g., Rybicki & Lightman 1979) and R(x) is given in Equation (23).

3.3.2. Synchrotron Self-Compton

Synchrotron radiation also serves as a source of seed photons for the SSC process, with corresponding spectrum

$$\begin{eqnarray}{f}_{{\epsilon }_{s}}^{\mathrm{SSC}} & = & \displaystyle \frac{9}{16}\displaystyle \frac{{(1+z)}^{2}{\sigma }_{{\rm{T}}}{\epsilon }_{s}^{{\prime} 2}}{\pi {\delta }_{{\rm{D}}}^{2}{c}^{2}{t}_{v}^{2}}{\displaystyle \int }_{0}^{\infty }\ d{\epsilon }_{* }^{{\prime} }\ \displaystyle \frac{{f}_{{\epsilon }_{* }}^{\mathrm{syn}}}{{\epsilon }_{* }^{{\prime} 3}}\ \\ & & \times {\displaystyle \int }_{{\gamma }_{1}^{{\prime} }}^{\infty }\ d{\gamma }^{{\prime} }\displaystyle \frac{{N}_{e}^{{\prime} }({\gamma }^{{\prime} })}{{\gamma }^{{\prime} 2}}{F}_{{\rm{C}}}\left(4{\gamma }^{{\prime} }{\epsilon }_{* }^{{\prime} },\displaystyle \frac{\epsilon }{{\gamma }^{{\prime} }}\right)\end{eqnarray} \tag{ 27 }$$

(e.g., Finke et al. 2008), where ${\epsilon }_{s}^{{\prime} }={\epsilon }_{s}(1+z)/{\delta }_{{\rm{D}}}$, ${\epsilon }_{* }^{{\prime} }={\epsilon }_{* }(1+z)/{\delta }_{{\rm{D}}}$, and

$$\begin{eqnarray}&&{\gamma }_{1}^{{\prime} }=\displaystyle \frac{1}{2}{\epsilon }_{s}^{{\prime} }\left(1+\sqrt{1+\displaystyle \frac{1}{{\epsilon }^{{\prime} }{\epsilon }_{s}^{{\prime} }}}\right).\end{eqnarray} \tag{ 28 }$$

The function ${F}_{{\rm{C}}}$ in Equation (27) is defined by

$$\begin{eqnarray}{F}_{{\rm{C}}}(p,q) & = & 2w\mathrm{ln}(w)+(1+2w)(1-w)\\ & & +\displaystyle \frac{{({pw})}^{2}}{2}\displaystyle \frac{1-w}{1+{pw}}\end{eqnarray} \tag{ 29 }$$

(Jones 1968; Blumenthal & Gould 1970), where

$$\begin{eqnarray}&&w=\displaystyle \frac{q}{p(1-q)}.\end{eqnarray} \tag{ 30 }$$

3.3.3. Compton-scattered External Radiation

The external Compton SED is calculated by assuming that each external photon source (whether BLR line or dust) is monochromatic and isotropic and homogeneous in the stationary frame, so that

$$\begin{eqnarray}{f}_{{\epsilon }_{s}}^{\mathrm{EC}} & = & \displaystyle \frac{3}{4}\displaystyle \frac{c{\sigma }_{{\rm{T}}}{\epsilon }_{s}^{2}}{4\pi {d}_{L}^{2}}\displaystyle \frac{{u}_{* }}{{\epsilon }_{* }^{2}}{\delta }_{{\rm{D}}}^{3}\\ & & \times {\displaystyle \int }_{{\gamma }_{1}}^{{\gamma }_{\max }}d\gamma \displaystyle \frac{{N}_{e}^{{\prime} }(\gamma /{\delta }_{{\rm{D}}})}{{\gamma }^{2}}{F}_{{\rm{C}}}\left(4\gamma {\epsilon }_{* },\displaystyle \frac{{\epsilon }_{s}}{\gamma }\right),\end{eqnarray} \tag{ 31 }$$

where

$$\begin{eqnarray}&&{\gamma }_{1}=\displaystyle \frac{1}{2}{\epsilon }_{s}\left(1+\sqrt{1+\displaystyle \frac{1}{\epsilon {\epsilon }_{s}}}\right)\end{eqnarray} \tag{ 32 }$$

(e.g., Georganopoulos et al. 2001; Dermer et al. 2009).

We apply the EC formalism to radiation from the disk, the dust torus, and the BLR. In this type of calculation, Dermer et al. (2014) assumed that Lyα will be the most prominent emission line in the BLR and therefore that all other emission lines can be ignored. We studied this assumption in detail by implementing one version of our model in which Lyα was the only emission from the BLR to undergo EC, in addition to dust (J = 2). We also implemented another version in which 25 additional lines from the Vanden Berk et al. (2001) template were included as sources of BLR seed photons (J = 27). This latter model comprises dust seed photons plus the 26 lines for which Finke (2016) provided radius, luminosity, and energy density relationships based on SDSS composite quasar spectra (Vanden Berk et al. 2001). By adopting the formalism of Finke (2016), we are able to avoid increasing the number of model free parameters when including additional BLR constituents.

Each source of EC seed photons is used to generate a corresponding cooling term in the electron transport equation (Equation (1)), and therefore the associated components of the SED are computed self-consistently in our model. We provide further details on the calculation method for each EC spectral component below.

3.3.4. External Compton of Disk Photons

The disk emits photons in the optical/UV due to the thermal radiation described in Section 3.1. In principle, any photons incident on the blob will experience Compton scattering by the blob electrons, and we estimate that effect for the disk photons by calculating the resulting $\nu {F}_{\nu }$ flux using

$$\begin{eqnarray}&&{f}_{\epsilon }^{\mathrm{EC}/\mathrm{Disk}}=\displaystyle \frac{{r}_{e}^{2}{\epsilon }_{s}^{2}{L}_{\mathrm{disk}}{\delta }_{{\rm{D}}}^{3}}{16\pi {r}_{\mathrm{blob}}^{2}{d}_{{\rm{L}}}^{2}{\epsilon }_{\mathrm{disk}}^{2}}{\int }_{{\bar{\gamma }}_{\mathrm{low}}}^{\infty }d\gamma \displaystyle \frac{{N}_{e}^{{\prime} }(\gamma /{\delta }_{{\rm{D}}})}{{\gamma }^{2}}\,\bar{{\rm{\Xi }}}\end{eqnarray} \tag{ 33 }$$

(e.g., Dermer et al. 2009), where we approximate the disk spectrum as monochromatic with an energy ${m}_{e}{c}^{2}{\epsilon }_{\mathrm{disk}}\approx 10\,\mathrm{eV}$. The function $\bar{{\rm{\Xi }}}$ in the integrand is defined by

$$\begin{eqnarray}&&\bar{{\rm{\Xi }}}\equiv 1-y+\displaystyle \frac{1}{1-y}-\displaystyle \frac{2y}{\bar{\epsilon }(1-y)}+{\left[\displaystyle \frac{y}{\bar{\epsilon }(1-y)}\right]}^{2},\end{eqnarray} \tag{ 34 }$$

where $y={\epsilon }_{s}/\gamma$ and $\bar{\epsilon }=\gamma {\epsilon }_{\mathrm{disk}}(1-{\mu }_{s})$. The lower bound of the integration in Equation (33) is

$$\begin{eqnarray}&&{\bar{\gamma }}_{\mathrm{low}}=\displaystyle \frac{{\epsilon }_{s}}{2}\left[1+\sqrt{1+\displaystyle \frac{2}{{\epsilon }_{\mathrm{disk}}{\epsilon }_{s}(1-{\mu }_{s})}}\,\right],\end{eqnarray} \tag{ 35 }$$

where ${\mu }_{s}=\cos \theta$. We assume $\theta =1/{\rm{\Gamma }}$ so that ${\rm{\Gamma }}={\delta }_{{\rm{D}}}$.

We implement this calculation, but find that the contribution of EC-scattering of disk photons is negligible in comparison to the rest of the SED components. Therefore, we neglect it in the remainder of the analysis.

3.3.5. External Compton of Dust Torus Photons

The dust torus emits an infrared distribution of blackbody photons, powered by the absorption of disk radiation. For the purpose of Compton scattering, we can approximate the dust emission spectrum as a monochromatic line with energy given by the peak of the blackbody. The corresponding dimensionless line energy for the dust temperature ${T}_{\mathrm{dust}}$ is

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{dust}}=5\times {10}^{-7}\left(\displaystyle \frac{{T}_{\mathrm{dust}}}{1000\ {\rm{K}}}\right),\end{eqnarray} \tag{ 36 }$$

and the associated energy density is given by

$$\begin{eqnarray}&&{u}_{\mathrm{dust}}=2.2\times {10}^{-5}\left(\displaystyle \frac{\xi }{0.1}\right){\left(\displaystyle \frac{{T}_{\mathrm{dust}}}{1000{\rm{K}}}\right)}^{5.2}\,\mathrm{erg}\ {\mathrm{cm}}^{-3}\ \end{eqnarray} \tag{ 37 }$$

(Nenkova et al. 2008; Sikora et al. 2009), where $\xi \equiv {L}^{\mathrm{dust}}/{L}^{\mathrm{disk}}\leqslant 1$ is the efficiency with which the dust reprocesses disk photons. In our calculations, we set $\xi =0.1$, which is an appropriate covering factor for the torus (Sikora et al. 2009). The dust temperature ${T}_{\mathrm{dust}}$ must be less than the dust sublimation temperature, and therefore ${T}_{\mathrm{dust}}\lesssim 2000$ K.

3.3.6. External Compton of BLR Emission Line Photons

The spectrum of the Compton-scattered BLR emission for each line depends crucially on the radius where the line photons are generated in the BLR (${r}_{\mathrm{line}}$) and on the location of the jet blob relative to the line source. Finke (2016) developed a simple way to estimate the location and emission properties of the various lines found in the composite SDSS quasar template (Vanden Berk et al. 2001) by using empirical relations derived from reverberation mapping. From Equation (16), we can find the disk luminosity for the photon energy $\epsilon$ corresponding to the wavelength$5100\,\mathring{\rm A}$, which we denote as ${L}_{5100\mathring{\rm A} }$. Reverberation mapping indicates that the distance of the Hβ source from the BH is then given by

$$\begin{eqnarray}&&{r}_{{\rm{H}}\beta }={10}^{16.94}{\left(\displaystyle \frac{{L}_{5100\mathring{\rm A} }}{{10}^{44}\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{0.533}\qquad \mathrm{cm}\end{eqnarray} \tag{ 38 }$$

(Bentz et al. 2013), neglecting uncertainties. The distance from the BH of the peak emission for a particular emission line, ${r}_{\mathrm{line}}$, is related to the distance of the Hβ line source from the BH by

$$\begin{eqnarray}&&{r}_{\mathrm{line}}={r}_{{\rm{H}}\beta }\,{\rho }_{\mathrm{line}},\end{eqnarray} \tag{ 39 }$$

where the quantity ${\rho }_{\mathrm{line}}$ is provided in Table 5 of Finke (2016) for all of the broad lines in the composite SDSS quasar spectrum of Vanden Berk et al. (2001). This table is predicated on the assumption that the ratios ${\rho }_{\mathrm{line}}$ are the same for all quasars.

Next, we note that the specific luminosity for the Hβ line is given by

$$\begin{eqnarray}&&{L}_{{\rm{H}}\beta }=1.425\times {10}^{42}\,\mathrm{erg}\ \,\ {{\rm{s}}}^{-1}{\left(\displaystyle \frac{{L}_{5100\mathring{\rm A} }}{{10}^{44}\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{1/0.8826}\ \end{eqnarray} \tag{ 40 }$$

(Greene & Ho 2005). The luminosity of a particular emission line, ${L}_{\mathrm{line}}$, is related to the luminosity of the Hβ line by

$$\begin{eqnarray}&&{L}_{\mathrm{line}}={L}_{{\rm{H}}\beta }\,{{\ell }}_{\mathrm{line}},\end{eqnarray} \tag{ 41 }$$

where the quantities ${{\ell }}_{\mathrm{line}}$ are tabulated in Table 5 of Finke (2016).

The energy density of a particular emission line can be approximated as

$$\begin{eqnarray}&&{u}_{\mathrm{line}}=\displaystyle \frac{{u}_{\mathrm{line},0}}{1+{({r}_{\mathrm{blob}}/{r}_{\mathrm{line}})}^{\beta }},\end{eqnarray} \tag{ 42 }$$

where

$$\begin{eqnarray}&&{u}_{\mathrm{line},0}=\displaystyle \frac{{L}_{\mathrm{line}}}{4\pi {{cr}}_{\mathrm{line}}^{2}}\end{eqnarray} \tag{ 43 }$$

(e.g., Hayashida et al. 2012; Dermer et al. 2014). If ${r}_{\mathrm{blob}}\ll {r}_{\mathrm{line}}$, then ${u}_{\mathrm{line}}\to {u}_{\mathrm{line},0}$. Detailed geometric calculations by Finke (2016) indicate that $\beta \approx 7.7$, which we use in all our calculations.

We use the method above to calculate the energy densities ${u}_{\mathrm{ph}}^{(j)}$ for all of the BLR seed photons and also for the dust seed photons. We use the energy densities to compute ${b}_{{\rm{C}}}$, the coefficient of the Compton-cooling term, and to calculate the observed $\nu {F}_{\nu }$ spectrum using Equation (31). An example of our full SED calculation with all our spectral components included is provided in Figure 3, which is further discussed in Section 5.

Zoom In Zoom Out Reset image size

We compute the stationary-frame jet powers resulting from our models to verify that there is enough energy available in the accretion flow plus the rotating supermassive BH to power the jet and evaluate whether or not it is in equipartition. The Poynting flux (magnetic field) jet power is computed using

$$\begin{eqnarray}&&{P}_{B}=2\pi {R}_{b}^{{\prime} 2}\beta c{{\rm{\Gamma }}}^{2}{u}_{B},\end{eqnarray} \tag{ 44 }$$

and the electron jet power is likewise given by

$$\begin{eqnarray}&&{P}_{e}=2\pi {R}_{b}^{{\prime} 2}\beta c{{\rm{\Gamma }}}^{2}{u}_{e},\end{eqnarray} \tag{ 45 }$$

where

$$\begin{eqnarray}&&{u}_{e}=\displaystyle \frac{{m}_{e}{c}^{2}}{{V}_{b}^{{\prime} }}\int d{\gamma }^{{\prime} }{\gamma }^{{\prime} }{N}_{e}^{{\prime} }({\gamma }^{{\prime} })\end{eqnarray} \tag{ 46 }$$

and ${V}_{b}^{{\prime} }=4\pi {R}_{b}^{{\prime} 3}/3$. The jet powers we obtain are reported in Table 4.

Hayashida et al. (2012) modeled the SEDs for 3C 279 for the same epochs that we analyze here. These authors utilized two models: a propagating emission region model and a jet precession scenario, which included two synchrotron components. They chose not to reproduce the X-rays with their models for most of the epochs, since the polarization and variability differences between the X-rays and the γ-rays suggested that they were produced by different electron populations. They usually assumed that the blob was outside the BLR. In most of our models, the X-rays come from SSC, the exceptions being the hard X-rays in Epochs A and C, which have a small contribution from EC. The SSC can vary with a changing ${R}_{b}^{{\prime} }$ independent of the optical (synchrotron) and $\gamma$-rays (EC) in our model.

Dermer et al. (2014) modeled the full SEDs, including the X-rays, but they made a key assumption that differentiates their model from ours: they assumed equipartition among the energy densities in the emitting region. They chose to reproduce the X-ray emission along with the rest of the IR through $\gamma$-ray SEDs. However, they had difficulty explaining the $\gamma$-ray data at $\gtrsim 5\,\mathrm{GeV}$. Our J = 2 models are most similar to the Dermer et al. (2014) models, since these only have EC/Lyα and EC/dust. For these J = 2 models, we still reproduce the data better than Dermer et al. (2014).

It is important to note that for a given electron distribution, our spectral calculations are identical to those of Dermer et al. (2014). Hence, the primary difference between our model and theirs is in how the electron distribution is determined. In our model, the electron distribution is determined self-consistently via the solution to the electron transport equation (Section 2.1) rather than using a simple ad hoc electron distribution. Furthermore, we do not assume equipartition. The assumption of equipartition gives Dermer et al. (2014) another constraint and less flexibility, which could account for why our models reproduce the data better. Our model electron distributions more closely resemble power laws with exponential cutoffs, due to the inclusion of the full Compton cross-section and the inefficient Compton cooling in the extreme KN regime. These effects allow the electrons to be accelerated to higher energies than they would be able to achieve in a Thomson model (e.g., Dermer & Atoyan 2002; Moderski et al. 2005; see Figure 2). Hence, the electrons in our model radiate higher energy $\gamma$-rays, in better agreement with the data.

Asano & Hayashida (2015) examined Epoch D with a time-dependent stochastic acceleration model in which the electrons are accelerated from some distance ${R}_{0}=0.023$ pc out to 2R₀ and then radiate. In this picture, $2{R}_{0}=1.4\times {10}^{17}$ cm, which is very similar to the blob distance ${r}_{\mathrm{blob}}$ that we find for Epoch D in our J = 2 model. This is consistent with their use of the same EC model as Hayashida et al. (2012), including external radiation from the BLR (as a single source) and dust. Their magnetic field B varies with the blob distance, but when acceleration turns off and radiation turns on, $B=3.5\,$ G, which is more than six times what we found in Epoch D. Their Doppler factor ${\delta }_{{\rm{D}}}=15$ is somewhat smaller, but comparable to our value of ${\delta }_{{\rm{D}}}=20$. Their energy densities are appreciably larger than ours (see Table 4), with ${u}_{B}=0.48$ erg cm⁻³and ${u}_{\mathrm{ext}}$ not provided directly but stated to be much larger. In contrast, our ${u}_{B}=1.2\times {10}^{-2}$ (erg cm⁻³) for J = 2 in Epoch D, and ${u}_{\mathrm{ext}}$ is almost 100 times smaller. So, their jet power must be very strongly field dominated even though their particle injection rate ${\dot{N}}_{e}=7.8\times {10}^{49}$ s⁻¹ of ${\gamma }_{e}=10$ electrons is $2\times {10}^{14}$ s⁻¹ times higher than our injection of ${\gamma }_{\mathrm{inj}}=1.01$ electrons.

Dermer et al. (2014) assumed equipartition between the energy densities of the particles and the magnetic field in order to constrain their SED modeling of 3C 279. Equipartition results in jet powers close to the minimum jet power. We do not assume equipartition in our modeling. We compute the equipartition parameter ${\zeta }_{e}={u}_{e}/{u}_{B}={P}_{e}/{P}_{B}$, following Dermer et al. (2014), and the results are reported in Table 4. Our models are either electron dominated or near equipartition. Since some of these results are out of equipartition, and hence far from the minimum power condition, we also compute the total jet power (${P}_{\mathrm{tot}}={P}_{e}+{P}_{B}$) as a fraction of accretion power (${P}_{\mathrm{acc}}={L}_{\mathrm{disk}}/0.4$, where the 0.4 is the efficiency for a maximally rotating BH), and these results are also reported in Table 4. In all but two of our models, we find that ${P}_{\mathrm{tot}}/{P}_{\mathrm{acc}}\lt 1$, indicating that accretion is sufficient to explain the power in the jet. General relativistic magnetohydrodynamic simulations carried out by Tchekhovskoy et al. (2011) indicate that ${L}_{\mathrm{tot}}/{P}_{\mathrm{acc}}\sim$ a few is possible with magnetically arrested accretion, which helps to support our high values. In the case of a spinning BH, significant energy can be drawn out of the ergosphere via the Blandford–Znajek process, reducing the need to power the blazar jet using accretion alone.

One of the strengths of utilizing a particle transport equation such as the one employed here (Equation (1)) is that the formalism automatically conserves energy. This provides us with a unique means for gaining further insight into how the energy is distributed among the various available channels associated with the acceleration, loss, radiative, and escape processes experienced by the electrons in the blob. We can use this information to explore the power distribution for each of the epochs of the 3C 279 data we analyze here and to determine what patterns may exist in the data that could provide significant clues about the relevant underlying physical processes. Furthermore, calculating the power in each channel also allows us to confirm that energy is conserved in the model, which lends further support to our conclusions. It is also interesting to confirm that the rate of electron escape from the blob is equal to the injection rate, which is non-trivial in our application since the escape timescale is energy dependent (see Equation (10)).

The steady-state transport equation with the full Compton cross-section implemented was solved using a numerical procedure, as described in Section 2.3 and the Appendix. The rate of escape of electrons from the blob, measured in the blob frame, is computed using

$$\begin{eqnarray}&&{\dot{N}}_{\mathrm{esc}}^{{\prime} }={\int }_{{\gamma }_{\min }^{{\prime} }}^{{\gamma }_{\max }^{{\prime} }}\displaystyle \frac{{D}_{0}{\gamma }^{{\prime} }}{\tau }{N}_{e}^{{\prime} }({\gamma }^{{\prime} })d{\gamma }^{{\prime} },\end{eqnarray} \tag{ 47 }$$

where the primes in this section emphasize that the powers, injection rates, and escape rates considered here are calculated in the co-moving (blob) frame, whereas the jet powers shown in Table 4 are given in the stationary (BH) frame. The integral in Equation (47) is computed using the full numerical solution for the electron distribution ${N}_{e}^{{\prime} }(\gamma ^{\prime} )$. In Table 5, we compare the results obtained for ${\dot{N}}_{\mathrm{esc}}^{{\prime} }$ with the electron injection rate, ${\dot{N}}_{\mathrm{inj}}^{{\prime} }$, and we note that for each epoch analyzed, we obtain identical values for these two quantities. This confirms that the electron "number budget" is correctly handled in our formalism. The particle injection rate is related to the injection luminosity ${L}_{{\rm{e}},\mathrm{inj}}^{{\prime} }$ via Equation (12).

Table 5.  Electron Energy Budget in the Co-moving Blob Frame

Parameter	Epoch (Model)
(Units)	A (J = 2)	A (J = 27)	B (J = 2)	B (J = 27)	C (J = 2)	C (J = 27)	D (J = 2)	D (J = 27)
(s−1)
${\dot{N}}_{\mathrm{inj}}^{{\prime} }$	$\quad 8.90\times {10}^{35}$	$\quad 8.80\times {10}^{35}$	$\quad 2.10\times {10}^{35}$	$\quad 1.35\times {10}^{35}$	$\quad 8.22\times {10}^{35}$	$\quad 5.00\times {10}^{35}$	$\quad 4.00\times {10}^{35}$	$\quad 3.85\times {10}^{35}$
${\dot{N}}_{\mathrm{esc}}^{{\prime} }$	$-8.90\times {10}^{35}$	$-8.80\times {10}^{35}$	$-2.10\times {10}^{35}$	$-1.35\times {10}^{35}$	$-8.22\times {10}^{35}$	$-5.00\times {10}^{35}$	$-4.00\times {10}^{35}$	$-3.85\times {10}^{35}$
(erg s−1)
${P}_{\mathrm{inj}}^{{\prime} }$	$\quad 7.36\times {10}^{29}$	$\quad 7.28\times {10}^{29}$	$\quad 1.74\times {10}^{29}$	$\quad 1.12\times {10}^{29}$	$\quad 6.80\times {10}^{29}$	$\quad 4.14\times {10}^{29}$	$\quad 3.31\times {10}^{29}$	$\quad 3.18\times {10}^{29}$
${P}_{\mathrm{esc}}^{{\prime} }$	$-4.21\times {10}^{31}$	$-4.51\times {10}^{31}$	$-2.35\times {10}^{31}$	$-1.35\times {10}^{31}$	$-2.31\times {10}^{30}$	$-1.24\times {10}^{31}$	$-7.36\times {10}^{31}$	$-6.81\times {10}^{31}$
${P}_{\mathrm{sto}}^{{\prime} }$	$\quad 7.64\times {10}^{42}$	$\quad 7.95\times {10}^{42}$	$\quad 8.00\times {10}^{42}$	$\quad 1.67\times {10}^{42}$	$\quad 4.23\times {10}^{41}$	$\quad 8.28\times {10}^{41}$	$\quad 1.50\times {10}^{44}$	$\quad 1.54\times {10}^{44}$
${P}_{\mathrm{sh},\mathrm{ad}}^{{\prime} }$	$-7.26\times {10}^{42}$	$-7.55\times {10}^{42}$	$-7.20\times {10}^{42}$	$-1.46\times {10}^{42}$	$-4.28\times {10}^{41}$	$-8.08\times {10}^{41}$	$-1.39\times {10}^{44}$	$-1.42\times {10}^{44}$
${P}_{\mathrm{syn}}^{{\prime} }$	$-6.97\times {10}^{40}$	$-6.60\times {10}^{40}$	$-1.31\times {10}^{41}$	$-3.98\times {10}^{40}$	$-3.96\times {10}^{38}$	$-4.62\times {10}^{39}$	$-1.45\times {10}^{42}$	$-1.49\times {10}^{42}$
${P}_{\mathrm{EC}}^{{\prime} }$	$-3.14\times {10}^{41}$	$-3.33\times {10}^{41}$	$-6.70\times {10}^{41}$	$-1.69\times {10}^{41}$	$-1.31\times {10}^{39}$	$-1.74\times {10}^{40}$	$-9.85\times {10}^{42}$	$-1.00\times {10}^{43}$
${P}_{\mathrm{net}}^{{\prime} }$	$-1.43\times {10}^{39}$	$-1.62\times {10}^{39}$	$-5.82\times {10}^{38}$	$-1.19\times {10}^{38}$	$-6.99\times {10}^{39}$	$-1.31\times {10}^{39}$	$-1.57\times {10}^{40}$	$-1.55\times {10}^{40}$
${\delta }_{\mathrm{err}}$	$0.02 \%$	$0.02 \%$	$0.01 \%$	$0.01 \%$	$1.65 \%$	$0.16 \%$	$0.01 \%$	$0.01 \%$

Download table as:  ASCIITypeset image

In addition to tracking the particle number budget, we can also calculate the power produced by each individual process in the transport equation in order to verify the energy budget in our numerical routine balances properly and to compare the relative energy contributions of each process affecting the electron distribution in the blob. The power provided by incident particles is the same as the injection luminosity (${P}_{\mathrm{inj}}^{{\prime} }={L}_{{\rm{e}},\mathrm{inj}}^{{\prime} }$). The escape power is given by

$$\begin{eqnarray}&&{P}_{\mathrm{esc}}^{{\prime} }={m}_{e}{c}^{2}\displaystyle \frac{{D}_{0}}{\tau }{\int }_{{\gamma }_{\min }^{{\prime} }}^{{\gamma }_{\max }^{{\prime} }}{{\gamma }^{{\prime} }}^{2}{N}_{e}^{{\prime} }({\gamma }^{{\prime} })d{\gamma }^{{\prime} }.\end{eqnarray} \tag{ 48 }$$

The remaining components of the electron power budget are calculated by integrating the associated terms in the Fokker–Planck drift coefficient (Equation (7)) since the broadening coefficient does not contribute to the energy budget. The powers associated with stochastic acceleration, shock acceleration+adiabatic losses, synchrotron losses, and inverse-Compton (EC) losses are denoted by ${P}_{\mathrm{sto}}^{{\prime} }$, ${P}_{\mathrm{sh},\mathrm{ad}}^{{\prime} }$, ${P}_{\mathrm{syn}}^{{\prime} }$, and ${P}_{\mathrm{EC}}^{{\prime} }$, respectively. Each of these quantities can be computed by integrating the respective term in the drift coefficient using the template

$$\begin{eqnarray}&&{P}_{\mathrm{drift}}^{{\prime} }={m}_{e}{c}^{2}{D}_{0}{\int }_{{\gamma }_{\min }^{{\prime} }}^{{\gamma }_{\max }^{{\prime} }}{C}_{\mathrm{drift}}({\gamma }^{{\prime} }){N}_{e}^{{\prime} }({\gamma }^{{\prime} })d{\gamma }^{{\prime} },\end{eqnarray} \tag{ 49 }$$

where

$$\begin{eqnarray}&&{C}_{\mathrm{drift}}({\gamma }^{{\prime} })=\left\{\begin{array}{l}4{\gamma }^{{\prime} }\\ a{\gamma }^{{\prime} }\\ -{b}_{\mathrm{syn}}{{\gamma }^{{\prime} }}^{2}\\ -{b}_{{\rm{C}}}{{\gamma }^{{\prime} }}^{2}\end{array}\right.\end{eqnarray} \tag{ 50 }$$

in turn for the cases of stochastic, shock/adiabatic, synchrotron, and Compton losses, respectively. The results obtained for each of these power components (in the blob frame) for each model and epoch of 3C 279 treated here are reported in Table 5. We also compute the sum of the gain and loss terms, given by

$$\begin{eqnarray}&&{P}_{\mathrm{net}}^{{\prime} }={P}_{\mathrm{inj}}^{{\prime} }+{P}_{\mathrm{esc}}^{{\prime} }+{P}_{\mathrm{sto}}^{{\prime} }+{P}_{\mathrm{sh},\mathrm{ad}}^{{\prime} }+{P}_{\mathrm{syn}}^{{\prime} }+{P}_{\mathrm{EC}}^{{\prime} }.\end{eqnarray} \tag{ 51 }$$

Ideally, we should obtain ${P}_{\mathrm{net}}^{{\prime} }=0$ if energy is perfectly conserved in our model. However, due to numerical errors, this can never be the case. We therefore compute a relative error, ${\delta }_{\mathrm{err}}$, by dividing $| {P}_{\mathrm{net}}^{{\prime} }|$ by the sum of the stochastic acceleration rate, ${P}_{\mathrm{sto}}^{{\prime} }$, and the injection rate ${P}_{\mathrm{inj}}^{{\prime} }$, which we report in Table 5. We confirm that the relative error is quite small in most cases, which is sufficient to confirm that energy is properly conserved in our model applications to each epoch of the 3C 279 data treated here.

Comparing the various channel powers listed in Table 5, we are able to draw some interesting conclusions about the physical processes influencing the behavior of the electrons in the blob. The particle injection power is always more than 10 orders of magnitude lower than the rest of the positive power contributions because the rate of particle injection is very low. As a result, the steady-state electron energy distributions we obtain are insensitive to the precise energy (or distribution) of the injected particles. The lack of significant pre-acceleration of the electrons is consistent with the observed lack of high-energy radiation along the upstream path of the jet, which implies that the precursor jet is "cold." In this scenario, the energy of the electrons in the precursor jet is mostly in the form of directed (bulk) kinetic energy rather than stochastic energy. This situation reverses once the electrons are energized in the acceleration region. We therefore conclude that the acceleration zone is located coincident with, or at least near, the emission region, which also helps to motivate the one-zone model considered here.

The stochastic acceleration power and the combined shock acceleration and adiabatic expansion power are the largest contributors to the overall energy budget. Stochastic acceleration provides the bulk of the acceleration to the electrons in the blob, while adiabatic expansion is the largest energy-loss mechanism. This is true for all epochs of 3C 279 that we examined, as can be seen in Table 5. The shock acceleration and adiabatic expansion processes are combined into a single term in our model, represented by the constant a in the transport equation, because they have the same first-order dependence on the particle energy. In general, adiabatic losses are expected to lead to the value $a\sim -4$ (Lewis et al. 2016). Since the values for a reported in Table 3 are all close to −4, this suggests that shock acceleration (which always makes a positive contribution to a) is a weak contributor to the acceleration of the electrons. Hence, we conclude that stochastic acceleration is the main contributor to the acceleration of the electrons in the blob.

Comparison of our model predictions with the observational data leads to the conclusion that stochastic acceleration and adiabatic expansion are the dominant processes in the examined epochs of 3C 279 (due to the breadth of the individual spectral components), as discussed in Section 4 and depicted in Figure 5. The active region, where the particle acceleration and emission takes place, is located $\sim 0.1$ pc from the central engine. We note that the generation of high-energy emission at this distance is consistent with Fermi LAT observations (e.g., Madejski & Sikora 2016). The occurrence of strong adiabatic cooling $\sim 0.1$ pc from the BH would imply substantial local acceleration of the bulk flow, which is supported by the observations (e.g., Homan et al. 2015). The acceleration could be due to passage through a recollimation shock that creates a nozzle in the flow (e.g., Cohen et al. 2015). The simultaneous dominance of stochastic particle acceleration at the same location may reflect the action of wave excitation at the recollimation shock (Marscher 2014). Taken together, these concepts support the scenario explored in the present paper, in which adiabatic losses and stochastic acceleration both dominate at a distance of ∼0.1 pc during the examined epochs of 3C 279.

The energy losses to the electrons due to the EC process are always greater in magnitude (more negative) than those due to the synchrotron process for all epochs of 3C 279 considered here. This is expected, because the part of the radiation spectrum due to EC always displays a higher amplitude in the $\nu {F}_{\nu }$ spectrum than that due to the synchrotron component for all of the 3C 279 epochs we examined. The electron energy losses due to radiation processes (EC and synchrotron) are always less than the losses due to adiabatic expansion.

Diltz & Böttcher (2014) studied electron evolution and emission in blazars using a Fokker–Planck equation that included second-order Fermi acceleration processes. However, they did not include a term representing first-order Fermi acceleration. Furthermore, they treated the BLR as a blackbody at a particular radius, not taking into account the various broad lines at different radii.

Asano et al. (2014) developed a particle acceleration and emission model that is somewhat similar to ours, although they considered Kolmogorov and Bohm stochastic acceleration, in addition to the hard-sphere prescription that we implement here. Asano & Hayashida (2015) applied the time-dependent model of Asano et al. (2014) to study FSRQ flares, and they also used the steady-state model to reproduce Epoch D of 3C 279 from Hayashida et al. (2012). In the time-dependent model of Asano & Hayashida (2015), the electrons are injected into a blob and then stochastically accelerated until the blob location reaches a predetermined distance from the BH, at which point the acceleration is turned off and the radiative emission processes take over. They include synchrotron, SSC, EC/dust, and EC/BLR as their radiative processes, and they also include adiabatic energy losses.

It is interesting to highlight the differences between the treatment of the BLR seed photons for the EC process considered here versus previous treatments. We employ a radially stratified model in which 26 individual species of seed photons are emitted from separate, discrete radial shells. Hayashida et al. (2012) and Asano & Hayashida (2015) employ a similar discrete shell geometry, except that they consider only one "equivalent" line component, instead of utilizing an explicit multicomponent BLR emission model like ours. In their approach, the equivalent line is a conceptual average of the multiple seed photon lines known to exist in the environment surrounding the BH. Asano et al. (2014) employ an isotropic (non-stratified) shell geometry. In order to justify the geometrical assumption in their model, the blob must be located inside the BLR, but there is no such restriction in our model.

The electron acceleration in the model of Asano & Hayashida (2015) is purely stochastic, and there is no explicit treatment of shock acceleration. Instead, shock acceleration is implemented approximately by invoking a population of pre-accelerated electrons with a prescribed energy ${\gamma }_{\mathrm{inj}}=100$ that is much higher than ours (see also Asano & Mészáros 2011). We improve upon the model of Asano & Hayashida (2015) by including shock acceleration explicitly as well as an energy-dependent Bohm particle escape term.

Diltz & Böttcher (2014) numerically solved a time-dependent transport equation, including synchrotron and Compton energy losses, hard-sphere stochastic acceleration, and energy-independent particle escape. They focused on the injection of electrons with a power-law energy distribution, which was intended to model the pre-acceleration of the electrons due to encounters with shocks in the jet. They used both the BLR and dust torus emission as seed photon sources for EC, but they modeled the BLR as a single blackbody at a single radius, whereas we employ a detailed stratified model for the BLR. They also assumed that the EC seed photon distribution is isotropic, and consequently their physical interpretation is constrained to the case where the blob is located inside the BLR radius. Conversely, our model accommodates blob locations both inside and outside the BLR. The Diltz & Böttcher (2014) model fits their sample data set well, and they are able to conclude that a decrease in the overall acceleration timescale is related to a correlation between the optical and $\gamma$-ray bands, which was observed during several 3C 279 flares.

Chen et al. (2014) used a multizone inhomogeneous and time-dependent model to investigate FSRQ emission regions. They implement continuous stochastic acceleration, radiative cooling, particle escape, and the injection of a pre-accelerated electron distribution. Thus, shock acceleration was not included in their transport equation. Their treatment of the EC process assumes that the incident radiation is generated in the dust torus only, and hence they neglect the BLR contribution. Their results suggest that different types of flares are associated with different magnetic field orientations. They conclude that SSC models are not preferred for FSRQs, in agreement with our findings.

We have developed a novel calculation of the electron distribution describing the electrons in a blob propagating through a blazar jet. The model is based on a steady-state electron transport equation. After solving the transport equation, we used the resulting electron distributions to self-consistently calculate the multiwavelength spectrum, using in turn either 2 or 27 sources of seed photons for the EC process. We applied our model to interpret and approximately fit the 2008–2009 spectral data for 3C 279, reported by Hayashida et al. (2012) and further analyzed by Dermer et al. (2014).

Our inclusion of two acceleration mechanisms in the electron transport equation distinguishes our work from previous models because the shape of the electron distribution is allowed to vary as different processes become dominant. Hence, the solutions we obtain for the electron distribution improve upon the previous ad hoc distributions, and they also improve upon those obtained by solving a transport equation that only includes one form of acceleration. For example, a transport equation where only stochastic acceleration is present will always produce a log-parabola electron distribution (Tramacere et al. 2011). Models that do not include shock acceleration explicitly usually assume that the injected electrons are pre-accelerated (e.g., Diltz & Böttcher 2014; Tramacere et al. 2011; Chen et al. 2014; Asano & Hayashida 2015). In contrast, we do not assume particle pre-acceleration and instead treat all acceleration processes explicitly using the transport equation. Hence, our model is a more fundamental, first-principles representation of the physical processes occurring in the blob. Our primary findings are listed below:

• 1.  
  A one-zone leptonic model is sufficient to match the data in all epochs examined.
• 2.  
  Stochastic acceleration dominates over shock acceleration in all of the calculations performed here.
• 3.  
  The distance of the blob from the BH is dependent on the structure and composition of the BLR. Specifically, the blob is calculated to be farther from the BH when the BLR is stratified.
• 4.  
  In our analysis of Epochs A–D from Hayashida et al. (2012), shock acceleration is present but overwhelmed by adiabatic losses.
• 5.  
  The utilization of a detailed, stratified model for the BLR leads to improved estimates for the model parameters.
• 6.  
  We calculate the matter and magnetic energy densities and conclude that the jet is not always in equipartition, but sometimes matter dominated.

We plan to use the new model developed here to analyze data for additional blazar γ-ray flares in future work.

When the full Compton cross-section is utilized in the transport equation (Equation (1)), an analytical solution is no longer possible. We therefore solve this equation using a bi-directional Runge–Kutta method. This approach requires two separate integrations of the transport equation, with one starting at the high-energy boundary, $\gamma ={\gamma }_{\max }$, and the other starting at the low-energy boundary, $\gamma ={\gamma }_{\min }$. In order to apply this method, we therefore need to specify appropriate boundary conditions at the two ends of the computational domain. We begin by writing Equation (1) in the flux-conservation form,

$$\begin{eqnarray}&&\displaystyle \frac{1}{{D}_{0}}\displaystyle \frac{\partial {N}_{e}}{\partial t}=0=-\displaystyle \frac{1}{{D}_{0}}\displaystyle \frac{\partial \dot{{ \mathcal N }}}{\partial \gamma }-\displaystyle \frac{{N}_{e}\gamma }{\tau }+\displaystyle \frac{{\dot{N}}_{\mathrm{inj}}\delta (\gamma -{\gamma }_{\mathrm{inj}})}{{D}_{0}},\end{eqnarray} \tag{ 52 }$$

where the electron transport rate in the energy space is defined by

$$\begin{eqnarray}&&\dot{{ \mathcal N }}\equiv {D}_{0}\left\{-{\gamma }^{2}\displaystyle \frac{\partial {N}_{e}}{\partial \gamma }+\phantom{\rule{}{3.0ex}}[2\gamma +a\gamma \right.\\ &&\,\left.\left.-{b}_{\mathrm{syn}}{\gamma }^{2}-{\gamma }^{2}\displaystyle \sum _{j=1}^{J}{b}_{{\rm{C}}}^{(j)}H(\gamma {\epsilon }_{\mathrm{ph}}^{(j)})\right]{N}_{e}\right\}.\end{eqnarray} \tag{ 53 }$$

At very low energies ($\gamma \to 0$), the electron transport rate, $\dot{{ \mathcal N }}$, must vanish, since no particles can be transported into the negative energy domain. Likewise, at sufficiently high energies, we expect that $\dot{{ \mathcal N }}\to 0$ because the electrons cannot be accelerated to arbitrarily high energies.

The maximum energy achieved by the electrons will be comparable to the exponential turnover energy in the electron distribution, ${N}_{e}(\gamma )$, which corresponds to a balance between acceleration and losses. This energy, which we refer to as the "equilibrium energy," ${\gamma }_{\mathrm{eq}}$, can be estimated by setting the Fokker–Planck drift coefficient (Equation (7)) equal to zero. At sufficiently high electron energies, the Compton cross-section declines exponentially due to the KN effect, and therefore synchrotron losses will dominate. The equilibrium energy is therefore given by

$$\begin{eqnarray}&&{\gamma }_{\mathrm{eq}}=\displaystyle \frac{4\,+\,a}{{b}_{\mathrm{syn}}},\end{eqnarray} \tag{ 54 }$$

Lewis et al. (2016). In general, we find that ${\gamma }_{\mathrm{eq}}$ exceeds the maximum Lorentz factor in our computational domain, ${\gamma }_{\max }$, and therefore we are justified in applying the zero-flux boundary condition at both the high- and low-energy boundaries in our Runge–Kutta code. Our zero-flux boundary condition can therefore be written as

$$\begin{eqnarray}&&\dot{{ \mathcal N }}(\gamma )=0,\qquad \left\{\begin{array}{l}\gamma ={\gamma }_{\min },\\ \gamma ={\gamma }_{\max }.\end{array}\right.\end{eqnarray} \tag{ 55 }$$

Combining Equations (53) and (55) yields an expression for the derivative ${{dN}}_{e}/d\gamma$, which can be written as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {N}_{e}}{\partial \gamma }=\left[\displaystyle \frac{2+a}{\gamma }-{b}_{\mathrm{syn}}-\displaystyle \sum _{j=1}^{J}{b}_{{\rm{C}}}^{(j)}H(\gamma {\epsilon }_{\mathrm{ph}}^{(j)})\right]{N}_{e}(\gamma ).\end{eqnarray} \tag{ 56 }$$

Equation (56) is applied at $\gamma ={\gamma }_{\min }$ and also at $\gamma ={\gamma }_{\max }$ to calculate the starting derivative once a boundary value for N_e has been specified, as discussed below.

One of the challenges of the Runge–Kutta implementation is that the equation is mathematically stiff, and it therefore requires extremely high grid resolution in order to obtain sufficient accuracy using standard methods. We therefore develop a novel alternative approach here that is based on ensuring a proper normalization for the global electron energy distribution, ${N}_{e}(\gamma )$.

The general solution for the electron energy distribution in our application is obtained by numerically integrating the homogeneous version of Equation (52) obtained when $\gamma \ne {\gamma }_{\mathrm{inj}}$, so that the source term is not active. We use a bi-directional technique to integrate Equation (52) away from the high-energy boundary at $\gamma ={\gamma }_{\max }$ and likewise away from the low-energy boundary at $\gamma ={\gamma }_{\min }$. The two fundamental solutions for N_e in the regions $\gamma \leqslant {\gamma }_{\mathrm{inj}}$ and $\gamma \geqslant {\gamma }_{\mathrm{inj}}$ are denoted by ${f}_{\mathrm{inc}}$ and ${f}_{\mathrm{dec}}$, where the subscripts refer to the "incrementing" and "decrementing" solutions, respectively. The derivatives at the two respective boundaries are each set using Equation (56). We set ${f}_{\mathrm{inj}}=1$ at $\gamma ={\gamma }_{\min }$ and ${f}_{\mathrm{dec}}=1$ at $\gamma ={\gamma }_{\max }$ without loss of generality, since the global solution for ${N}_{e}(\gamma )$ will be renormalized later in the process.

Once the fundamental solutions ${f}_{\mathrm{inc}}$ and ${f}_{\mathrm{dec}}$ have been obtained, the normalized global solution for the electron number distribution is given by

$$\begin{eqnarray}&&{N}_{e}(\gamma )=\left\{\begin{array}{l}{{Af}}_{\mathrm{inc}}(\gamma ),\quad \gamma \leqslant {\gamma }_{\mathrm{inj}},\\ {{Bf}}_{\mathrm{dec}}(\gamma ),\quad \gamma \geqslant {\gamma }_{\mathrm{inj}},\end{array}\right.\end{eqnarray} \tag{ 57 }$$

where A and B are normalization coefficients that can be computed by ensuring a steady-state balance between electron injection and escape. The rate at which electrons escape from the blob is given by

$$\begin{eqnarray}&&{\dot{N}}_{\mathrm{esc}}={\int }_{0}^{\infty }\displaystyle \frac{{N}_{e}(\gamma )}{{t}_{\mathrm{esc}}(\gamma )}d\gamma ,\end{eqnarray} \tag{ 58 }$$

and the escape timescale, ${t}_{\mathrm{esc}}$, is related to the dimensionless parameter τ via (Lewis et al. 2016)

$$\begin{eqnarray}&&\displaystyle \frac{1}{{t}_{\mathrm{esc}}(\gamma )}=\displaystyle \frac{\gamma {D}_{0}}{\tau }.\end{eqnarray} \tag{ 59 }$$

In a steady state, the rate of escape of electrons from the blob, ${\dot{N}}_{\mathrm{esc}}$, must equal the injection rate, ${\dot{N}}_{\mathrm{inj}}$, computed using Equation (58). Hence, we require that

$$\begin{eqnarray}&&{\dot{N}}_{\mathrm{esc}}={\dot{N}}_{\mathrm{inj}}.\end{eqnarray} \tag{ 60 }$$

Combining relations, we find that

$$\begin{eqnarray}&&{\dot{N}}_{0}=A{\int }_{0}^{{\gamma }_{\mathrm{inj}}}\displaystyle \frac{{f}_{\mathrm{inc}}(\gamma )}{{t}_{\mathrm{esc}}(\gamma )}d\gamma +B{\int }_{{\gamma }_{\mathrm{inj}}}^{\infty }\displaystyle \frac{{f}_{\mathrm{dec}}(\gamma )}{{t}_{\mathrm{esc}}(\gamma )}d\gamma .\end{eqnarray} \tag{ 61 }$$

The electron distribution must be continuous at the injection energy, therefore

$$\begin{eqnarray}&&{{Af}}_{\mathrm{inc}}({\gamma }_{\mathrm{inj}})={{Bf}}_{\mathrm{dec}}({\gamma }_{\mathrm{inj}}).\end{eqnarray} \tag{ 62 }$$

Solving for B yields

$$\begin{eqnarray}&&B\,=\,A\,\displaystyle \frac{{f}_{\mathrm{inc}}({\gamma }_{\mathrm{inj}})}{{f}_{\mathrm{dec}}({\gamma }_{\mathrm{inj}})}.\end{eqnarray} \tag{ 63 }$$

We can eliminate B between Equations (61) and (63) and solve for A to obtain

$$\begin{eqnarray}&&A={\dot{N}}_{0}{\left[{\int }_{0}^{{\gamma }_{\mathrm{inj}}}\displaystyle \frac{{f}_{\mathrm{inc}}(\gamma )}{{t}_{\mathrm{esc}}(\gamma )}d\gamma +\displaystyle \frac{{f}_{\mathrm{inc}}({\gamma }_{\mathrm{inj}})}{{f}_{\mathrm{dec}}({\gamma }_{\mathrm{inj}})}{\int }_{{\gamma }_{\mathrm{inj}}}^{\infty }\displaystyle \frac{{f}_{\mathrm{dec}}(\gamma )}{{t}_{\mathrm{esc}}(\gamma )}d\gamma \right]}^{-1}.\end{eqnarray} \tag{ 64 }$$

The fully normalized global solution for the electron distribution, valid over the entire energy range, ${\gamma }_{\min }\leqslant \gamma \leqslant {\gamma }_{\max }$, can now be written as

$$\begin{eqnarray}{N}_{e}(\gamma )={{Af}}_{\mathrm{inc}}({\gamma }_{\mathrm{inj}})\left\{\begin{array}{ll}{f}_{\mathrm{inc}}(\gamma )/{f}_{\mathrm{inc}}({\gamma }_{\mathrm{inj}}), & \gamma \leqslant {\gamma }_{\mathrm{inj}},\\ {f}_{\mathrm{dec}}(\gamma )/{f}_{\mathrm{dec}}({\gamma }_{\mathrm{inj}}), & \gamma \geqslant {\gamma }_{\mathrm{inj}},\end{array}\right.\end{eqnarray} \tag{ 65 }$$

We employ a trapezoid numerical integration scheme to compute the constant A in Equation (25) based on the fundamental solutions ${f}_{\mathrm{inc}}$ and ${f}_{\mathrm{dec}}$.

The normalization technique developed here is theoretically equivalent to the usual Wronskian method, but the latter method requires a very accurate numerical computation of the Wronskian of the fundamental solutions ${f}_{\mathrm{inc}}$ and ${f}_{\mathrm{dec}}$. Our alternative normalization technique is much more accurate and efficient when dealing with stiff ODEs such as our transport equation.