A multi-zone model for simulating the high energy variability of TeV blazars
Philip B. Graff, Markos Georganopoulos, Eric S. Perlman, Demosthenes Kazanas
aa r X i v : . [ a s t r o - ph ] A ug A multi-zone model for simulating the high energy variability ofTeV blazars
Philip B. Graff , Markos Georganopoulos , , Eric S. Perlman , Demosthenes Kazanas ABSTRACT
We present a time-dependent multi-zone code for simulating the variabilityof Synchrotron-Self Compton (SSC) sources. The code adopts a multi-zone pipegeometry for the emission region, appropriate for simulating emission from astanding or propagating shock in a collimated jet. Variations in the injection ofrelativistic electrons in the inlet propagate along the length of the pipe coolingradiatively. Our code for the first time takes into account the non-local, time-retarded nature of synchrotron self-Compton (SSC) losses that are thought tobe dominant in TeV blazars. The observed synchrotron and SSC emission is fol-lowed self-consistently taking into account light travel time delays. At any giventime, the emitting portion of the pipe depends on the frequency and the natureof the variation followed. Our simulation employs only one additional physicalparameter relative to one-zone models, that of the pipe length and is computa-tionally very efficient, using simplified expressions for the SSC processes. Thecode will be useful for observers modeling GLAST, TeV, and X-ray observationsof SSC blazars.
Subject headings: galaxies: active — quasars: general — radiation mechanisms:nonthermal — X-rays: galaxies
1. Introduction
In blazars, radio loud active galaxies with their relativistic jets pointing close to our lineof sight (Blandford 1978), the observed radiation is dominated by relativistically beamed Department of Physics, Joint Center for Astrophysics, University of Maryland, Baltimore County, 1000Hilltop Circle, Baltimore, MD 21250, USA NASA Goddard Space Flight Center, Code 663, Greenbelt, MD 20771, USA Department of Physics and Space Sciences, Florida Institute of Technology, 150 West University Boule-vard, Melbourne, FL 32901, USA ∼ TeV energies, SSC is the dominantemission mechanism. Recent observational results (e.g. D’Arcangelo et al. 2007; Marscheret al. 2008), however, place the blazar emission site beyond the broad line region, lendingsupport to the possibility that even in powerful blazars the GeV emission process may bepure SSC. For a review of leptonic models, as well as hadronic models for blazar emission(e.g. Aharonian 2000) see B¨ottcher (2007).Due to the small angular size of the blazar emission region, it is not possible to spatiallyresolve the emitting region. Because of this, information about the structure of the emit-ting source can be obtained only through multiwavelength variability studies. Particularlytelling is the variability of the emission produced by the highest energy electrons, becausethese electrons lose energy very quickly and exist only close to the sites where they have beenproduced. The goal of multiwavelength variability campaigns, involving in many cases ob-servations from radio up to TeV energies, is to study the characteristics of blazar variability,such as correlations and/or time delays between different energies, spectral characteristics ofthe observed variability, and the amplitude of variability as a function of energy.Most notable amongst blazars are the so called TeV blazars for which the synchrotronemission peaks at X-ray energies and the SSC emission peaks at TeV energies, as theypresent the active galaxies producing the highest confirmed electron energies. Variationsof TeV blazars in these two bands can be extremely rapid (TeV doubling times as shortas a few min; Aharonian et al. 2007), suggesting highly relativistic sub-pc scale flows(Doppler factors δ ∼
50; e.g. Begelman, Fabian & Rees 2008) that decelerate substantially(Georganopoulos & Kazanas 2003; Ghisellini, Tavecchio & Chiaberge 2005) to match themuch slower speeds required by VLBI observations (Piner & Edwards 2004; Piner, Pant, Edwards2008). The TeV and X-ray variations are usually well correlated (e.g. Fossati et al. 2008;Maraschi et al. 1999; Sambruna et al. 2000), as expected, because they present variationsby the same electron population. Usually, the lower energy emission within each of thesebands peaks with small time delays relative to the higher energy emission (e.g. Fossati et al. 3 –2000), while the X-ray and TeV spectra become harder with increasing flux (e.g. Takahashiet al. 1996). In certain cases, however, the X-ray and TeV variability do not seem to becorrelated in a simple way (e.g. Aharonian et al. 2005). An intriguing variability patternis that of the so-called ‘orphan’ flares, rare TeV flares that are not accompanied by X-rayflares (e.g. Krawczynski et al. 2004, B la˙zejowski et. al. 2005). While the correlated X-ray- TeV flares can be understood through an increase of the high energy emitting electrons,orphan TeV flares defy such a straightforward explanation.Models of blazar emission to date have, for the most part, been in some form of homo-geneous one zone models (e.g. Mastichiadis & Kirk 1997; Krawczynski, Coppi, & Aharonian2002). Such models, although appropriate for modeling the steady-state emission of a source,cannot simulate variability faster than the zone light crossing time. The basic limitation ofone zone models stems from the fact that the high energy variability of both the synchrotronand SSC components is produced by high energy electrons with cooling times shorter thanthe light crossing time. Even if we assume that a disturbance in the radiating plasma (e.g.a higher density) instantaneously propagates across the zone, the received radiation wouldbe smeared out for timescales shorter than the light crossing time, due to light travel timedelays from different parts of the source ( § i ) the SSC emission inside any one of their single zones uses as seed photons only thesynchrotron photons produced in that zone and ( ii ) the SSC energy losses in every zone arecaused only by the synchrotron photons produced in that zone. This simplified approach is agood approximation for following the energetics of the electrons if the source is synchrotron 5 –dominated (because the SSC losses, although inappropriately calculated, are negligible), butdoes not produce realistic SSC light curves, because it does not calculate the emission dueto upscattering synchrotron photons produced in other parts of the source in retarded times.A significant improvement was introduced by Sokolov, Marscher, & McHardy (2004)who incorporated in the calculation of the SSC emission from a given location in an inhomo-geneous source the synchrotron photons produced throughout the source in retarded times.This produces accurate SSC light curves, provided that the SSC losses that were still treatedas a local process are negligible. In a follow up paper, Sokolov & Marscher (2005), consid-ered also external Compton photons from the broad line region and the molecular torus. Thechallenge for inhomogeneous multi zone models for sources such as the TeV emitting blazarsis the calculation of the non-local, time-retarded SSC losses induced by photons producedin other parts of the source.Here, we present such an inhomogeneous model that, for the first time, takes intoaccount the non-local, time-delayed source emission on the SSC losses. We assume that apower law of relativistic electrons is injected at the inlet of a pipe, and that the electrons flowdownstream and cool radiatively. Variations in the injected electron distribution propagatedownstream and manifest themselves as frequency dependent variability. This allows us tomodel high energy multiwavelength variability in a self-consistent manner. In § § § §
2. The One Zone Model
We consider a homogeneous spherical source of radius R permeated by a magnetic field B of energy density B / (8 π ). Energetic electrons are injected into the region at a rate q ( γ, t ), where γ is the electron Lorentz factor and t the injection time. These electrons loseenergy through synchrotron and inverse Compton radiation and eventually escape after acharacteristic time, t esc , of the order of the light crossing time. The implementation wedescribe is applicable to sources that are optically thin both to synchrotron emission in the 6 –frequency range under consideration and to γ -ray absorption due to pair-production.The kinetic equation that describes the time-evolution of the electron energy distribution(EED) n ( γ, t ) is ∂n ( γ, t ) ∂t + ∂∂γ [ ˙ γ n ( γ, t )] + n ( γ, t ) t esc = q ( γ, t ) . (1)Here, ˙ γ includes both the synchrotron losses ˙ γ s and the inverse Compton losses ˙ γ IC in theThomson regime ( ǫγ ≤ / γ = ˙ γ s + ˙ γ IC , ˙ γ s = 4 σ τ mc γ U B , ˙ γ IC = 4 σ τ mc γ Z min [ ǫ max , / (4 γ )] ǫ min U ( ǫ, t ) dǫ (2)where U ( ǫ, t ) is the photon field energy density, ǫ is the photon energy in units of theelectron rest energy m e c , and σ τ is the Thomson cross section. The photon field U ( ǫ, t )includes not only the synchrotron produced photons, but all the photons produced in thesource through IC scattering, including, therefore, all the higher order SSC emission.We calculate the synchrotron emission following Melrose (1980): L s ( ǫ s , t ) = 1 . √ q Bh Z γ max γ min z / e − z n ( γ, t ) dγ, z = (cid:18) (cid:19) / ǫ s /B ⋆ γ , (3)where q is the electron charge, B ⋆ = B/B crit , and B crit = ( m e c ) / ( e ~ ) is the critical magneticfield, where the electron cyclotron energy equals its rest mass and strong field considerationsbecome important (e.g. Harding & Lai 2006). We note, that, although a δ -function approachfor calculating the synchrotron emissivity would be faster, it would misinterpret the spectra incases of hard power - law injection q ∝ γ − p , p >
2, whose cooling is known to produce a pile-up at its high energy cutoff (e.g. Kardashev 1962). Also, a δ -function synchrotron emissivitywould not produce the f ν ∝ ν / spectrum at frequencies below the critical frequency of thelowest energy electrons. These lower energy photons can be important seed photons forproducing hard SSC TeV emission as Katarzy¨nski et al. (2006) point out.To obtain the SSC emission through a simple integration as in the synchrotron case, weemploy the δ -function approximation in which seed photons of energy ǫ are IC scattered byelectrons of Lorentz factor γ to energy ǫ IC = (4 / ǫ γ , as long as the scattering takes placein the Thomson regime ( γǫ < / dǫ beingIC scattered by electrons with Lorentz factors in the range dγ , then the emitted IC power is˙ γ IC m e c n ( γ, t ) dγ and it is spread over a final photon energy range dǫ IC = 4( ǫ ǫ IC / / dγ ,resulting in an IC specific luminosity per seed photon energy interval dL IC ( ǫ IC , t ) = m e c n ( γ, t ) ˙ γ IC dγdǫ IC δ ( γ − (3 ǫ IC / ǫ ) / )Θ(3 / − γǫ ) , (4) 7 –where in this context ˙ γ IC = (4 σ τ / mc ) γ U ( ǫ , t ) dǫ and Θ( x ) is the Heaviside step function.This is written as dL IC ( ǫ IC , t ) = 3 / σ τ cn ( γ, t ) U ( ǫ , t ) dǫ ǫ / IC ǫ / δ ( γ − (3 ǫ IC / ǫ ) / )Θ(3 / − γǫ ) . (5)Integrating over the available seed photon distribution we obtain L IC ( ǫ IC , t ) = 3 / σ τ cǫ / IC Z ǫ ,max ǫ ,min n ( γ, t ) U ( ǫ , t ) ǫ − / δ ( γ − (3 ǫ IC / ǫ ) / ) dǫ . (6)The range of final photon energies ǫ IC is (4 / ǫ seed,min γ min < ǫ IC < γ max . For ǫ IC withinthis range the limits of the above integration are: ǫ ,min = ǫ seed,min for ǫ IC ≤ ǫ seed,min γ max ǫ IC γ max for ǫ IC ≥ ǫ seed,min γ max , (7) ǫ ,max = ǫ IC γ min for ǫ IC ≤ γ min ǫ IC for ǫ IC ≥ γ min . (8)Following Chang & Cooper (1970) and Chiaberge & Ghisellini (1999), we discretisethe kinetic equation (1), using a grid of logarithmically spaced Lorentz factors, γ j , j =0 , , , ..., j max , and linearly spaced time indices, t i . The difference equation that describesthis system is n j,i +1 − n j,i ∆ t = − ˙ γ j +1 ,i +1 n j +1 ,i +1 − ˙ γ j,i +1 n j,i +1 ∆ γ + q j,i +1 − n j,i +1 t esc . (9)Note that this is an implicit scheme, in the sense that the calculation of n j,i +1 requiresknowledge not of only the previous timestep EED, but also n j +1 ,i +1 , the next higher γ gridpoint at the current time. It is due to the implicit nature of the numerical procedure thatthis scheme is stable for large timesteps. The difference equation can be written as a systemof tridiagonal equations n j,i +1 = a n j,i − b n j +1 ,i +1 + c q j,i +1 , (10) a = ∆ γ ∆ γ + ∆ t ∆ γ /t esc − ∆ t ˙ γ j,i +1 , b = a ∆ t ∆ γ ˙ γ j +1 ,i +1 , c = a ∆ t. (11) 8 –This can be easily computed if we use the initial condition n j, = 0 ∀ j (start with norelativistic electrons in the system) and the boundary condition n j max ,t = 0 ∀ t . The firstcondition implies that the initial photon field is also zero for all photon energies. Thesecond condition is satisfied if we set γ j max > γ max , because the electrons can only loseenergy, and there is no way to move to higher energies, populating the j max bin of the γ -grid. The simulation proceeds in the following manner: given n ( γ, t i ) and U ( ǫ , t i ) wefirst calculate n ( γ, t i +1 ). We then calculate the synchrotron luminosity, L S ( ǫ, t i +1 ), and theinverse Compton luminosity, L IC ( ǫ, t i +1 ). The specific photon energy density U ( ǫ , t i +1 ) forthe next timestep is obtained by adding these two luminosities and dividing by 4 πR m e c . By construction, in the one zone model variations in the injection propagate instanta-neously throughout the source, because no spatial coordinate enters the description of thesystem. If a power - law EED, q ∝ γ − p , γ min ≤ γ ≤ γ max , is injected in the source, radia-tive cooling and electron escape will result in a broken power law EED n ( γ ) in the source,steepening from an electron index p to p + 1, above γ b , the electron Lorentz factor for whichthe escape time equals the radiative loss time γ b ˙ γ = t esc ⇒ γ b = 3 m e c σ τ U t esc , (12)where U stands for the total photon and magnetic field energy density in the source. Itis these electrons, with γ > γ b that produce the high energy synchrotron and IC emission.Because the electron escape time is of the order of the light crossing time, t esc = kt lc = kR/c , k ∼ − few, electrons with Lorentz factor γ > γ lc = kγ b , have a cooling time shorter thanthe light crossing time. One therefore anticipates that even for an injection event lastingmuch less that t lc , the high energy variable emission produced by electrons with γ > γ lc willbe smeared out by light travel time effects and would appear to last for ∼ t lc , even thoughin each point in the source it lasts a shorter time ∼ t lc γ lc /γ .To demonstrate this, a flaring state was simulated by an increase by a factor of 5 inthe injection q ( γ, t ) that lasted t inj = t lc /
10. The system was allowed to reach a steadystate before the disturbance in the injection was introduced. The emitted luminosity as afunction of frequency was followed in time allowing us to produce light curves. By not takinginto account time delays, and wrongly assuming that at any given time the observer seesthe entire source as being at a single physical state, electrons with t cool < t lc produce highenergy synchrotron and SSC variations that last less than t lc (upper panel of Figure 1).To take time delays into account, one has to consider that if at a certain time t the 9 –observer receives photons from the nearest (front) part of the source, slices further awayfrom the observer will be seen as they were in retarded times t − r/c , where r is the distanceof the slice from the front part of the source. To treat this, at each time the luminosityof the system was recorded for a number of time steps covering a time equal to the lightcrossing time of the region. The luminosity observed at any time is thus the sum of theluminosity from each of these time steps, as each one represents the light emitted by a slicesequentially further back from the observer. The resulting light curves, plotted in the lowerpanel of Figure 1, show that when the size of the region, and thus the time taken by lightto travel across it, is accounted for, the observed variability of high energy electrons with t cool < t lc is spread out over the length of the light crossing time. This is similar to whatChiaberge & Ghisellini (1999) observed when they performed a similar test.Another serious problem stemming from the lack of spatial considerations in the one-zone model comes from the fact that the model by construction assumes that the photonsproduced in the source at a given time are instantaneously available as seed photons for ICscattering throughout the source. This unphysical assumption has serious implications onthe calculation of the SSC emissivity and on the calculation of the SSC losses, which in turnaffect the evolution of the EED in the source, and through this the entire spectrum and lightcurves. The first effect has been addressed by the inhomogeneous model of Sokolov et al.(2004) and Sokolov & Marscher (2005). We present now the first multi-zone simulation thatincorporates the issue of light travel time effects on the SSC losses.
3. The Multi Zone Model3.1. The flow geometry
The simplest and least computationally intensive deviation from a homogeneous modelthat can address the issues discussed above is one in which plasma is injected into a pipeof radius R and length L and cools radiatively as it flows downstream before it escapesafter traversing the pipe length. Physically, this resembles the situation of a standing orpropagating shock, as seen at the frame of the shock front. The plasma flow velocity u andthe magnetic field B are constant along the pipe, and the EED is assumed to have no lateralgradients along the cross section of the pipe. Relativistic plasma is injected at the base ofthe flow. The injection variation timescale in this geometry can be arbitrarily smaller than R/c without violating causality, because, in principle, a disturbance in the plasma flow canreach the entire cross-section of the inlet at a single instance. However, the discretizationprocedure we describe below limits the range of meaningful injection variability to timescalesgreater than the plasma flow time through a zone of the flow. Our goal is to calculate the 10 –EED in the frame of the pipe as a function of time and distance z from the inlet of the pipe,and through this calculate the emission received by an observer located at an angle θ to theaxis of the pipe. The pipe is broken down lengthwise and all cells are of length l , comparable to the piperadius R . The length of the pipe L = N l , where N is the number of cells. Each cell is thensimulated by a one zone model. The electron injection at the first cell is q ( γ, t ), similar tothat defined for the one zone model. In each time step, the electrons that are calculated toleave each cell, are injected into the next cell in line in the next timestep. The injection ofelectrons, therefore, in cell i at time t j is q i ( γ, t j ) = n i − ( γ, t j − ) /t esc and the kinetic equationfor the i -th cell is ∂n i ( γ, t j ) ∂t + ∂∂γ [ ˙ γ n i ( γ, t j )] + n i ( γ, t j ) t esc = n i − ( γ, t j − ) t esc (13)Because in a time ∼ t esc the electron content of a cell is transferred to the next cell, t esc is connected to the bulk flow velocity u through t esc = l/u . We also use t esc as the timestep of our simulation. This ensures that the actual distance a disturbance in the electrondistribution travels in a time step is equal to the bulk velocity times the time-step size, bytransferring in a time step the electron content of cell i to cell i + 1. The shortest variabilitytimescale that can be simulated by this configuration is the single cell escape time. Thislimits the highest energy electrons that can be followed accurately to that of Lorentz factor γ b = 3 m e cu/ σ τ U l , where U is the photon plus magnetic field energy density in the first cell(see equation 12), and through this the highest energies of synchrotron and IC variationsthat can be reproduced. The advantage of the pipe configuration relative to a homogeneousmodel of the same size is that the highest energy electron variability we can follow is notconnected to the length of the entire pipe, but to the length of a single zone, a quantitythat is N times shorter. This results in the pipe being able to track variations faster by N ,following electrons more energetic by N , and synchrotron and SSC fequencies higher by N ,relative to a homogeneous model of size L . An early version of this approach was presentedby Graff et al. (2007). To solve the kinetic equation for each cell, an expression for the photon energy densityresulting from all other cells by taking into account light travel time delays is required. In 11 –general, for a region S characterized by a time-dependent emission coefficient j ( r ′ , t ′ , ǫ ), thephoton energy density U ( r , t, ǫ ) is calculated by integrating j ( r ′ , t ′ , ǫ ) /c in retarded timesover the volume of the region S . Setting r = 0, for a point of interest in S , yields, withoutloss of generality, U ( r = 0 , t, ǫ ) = 1 c Z r ′ (Ω)0 j ( r ′ , t ′ = t − r ′ /c, ǫ ) dr ′ d Ω . (14)For our geometry we express this through the following approximation. Consider a cell i centered at z i being illuminated by a cell m centered at z m . The solid angle subtended at z i by the cell m is ∆Ω ≈ πR / ( z i − z m ) . The photon energy density ∆ U ( z i , z m , t, ǫ ) at ( z i , t)due to photons produced at z j at retarded time t ′ = t − | z i − z m | /c is:∆ U ( z i , z m , t, ǫ ) = j ( z m , t ′ = t − | z i − z m | /c, ǫ ) c πR ( z i − z m ) l. (15)Making use of the fact that the volume of each cell is V c = πR l , and that the luminosity L ( ǫ ) emitted from a cell is L ( ǫ ) = 4 πj ( ǫ ) V c , we obtain∆ U ( z i , z m , t, ǫ ) = L ( z m , t ′ = t − | z i − z m | /c, ǫ )4 πc ( z i − z m ) . (16)A summation over all cells in the pipe results to the total photon energy density U ( z i , t j , ǫ k )at cell i , time t j = jt esc , and energy ǫ k U ( z i , t j , ǫ k ) = L ( z i , t j , ǫ k )2 πcR ( R + l ) + N X m =0 ,m = i L ( z m , t ′ = jt esc − | i − m | l/c, ǫ k )4 πcl ( i − m ) , (17)where the first term is the photon energy density due to cell i itself. Note that a calculationof the photon energy density requires keeping record of a data-cube of the SED emitted byeach cell at each time for a number of time steps equal to the light crossing time of the entireregion. When calculating the observed luminosities, we must take into account the differentdistances that light must travel from each of the different cells in the pipe to the observer.This difference is a function of the angle θ formed between the pipe and the observer. Ifphotons emitted from the first cell are received by the observer at a given time, photons fromcell i that are received simultaneously were emitted earlier by i l cos θ/c . Beaming can be 12 –easily included in this model, by assuming that the entire pipe is moving with a relativisticvelocity along its axis. Then for a choice of bulk Lorentz factor Γ and orientation angle θ obs in the observer’s frame, one can calculate the Doppler factor δ and the angle θ in the frameof the pipe, and transform the arrival times by dividing by δ , the observed frequencies bymultiplying by δ , and the observed fluxes by multiplying by δ .
4. Results4.1. Comparison with analytical results for the synchrotron dominated case.
A simple analytical test can be performed in the case of a source in which the energylosses are dominated by synchrotron radiation. In this case, adopting a steady power lawelectron injection at the inlet and assuming an electron residence time kL/c in the source,the source integrated EED will reach after time t >> kL/c a steady-state. This steady-state, source-integrated electron distribution is a broken power law with an electron indexsteepening by one at γ b = (3 m e c ) / (4 σ τ U B kL ). This will produce a synchrotron spectrumwith a spectral break of 1 / ǫ b = B ⋆ γ b . For the synchrotron dominatedconfiguration presented in Figure 2, the numerical result is in good agreement with theanalytical both for the EED and the SED employing a δ -function synchrotron emissivity.Note that the SSC luminosity is much lower than the synchrotron one. Note also that whilethe analytical synchrotron emission stops at ǫ min = B ⋆ γ min , the numerical continues to lowerfrequencies, due to the f ν ∝ ν / lower energy tail of the synchrotron emissivity.To compare the variability of our code with analytic expectations, we initiate injection ofa power low EED at t = 0 at the inlet and follow the evolution of the system toward a steadystate. We expect that the EED will reach a steady state at or before a time ∼ kL/c = 2 L/c ,equal to the time it takes for an electron to transverse the length of the pipe and escape,or equivalently the time it takes to fill the pipe with electrons. The part of the EED with γ < γ b will reach a steady state at ∼ kL/c , because the electrons responsible for this emissiontransverse the entire pipe without cooling appreciably. At higher energies, γ > γ b , theelectron radiative cooling lifetime t = (3 m e c ) / (4 σ τ U B γ ) is shorter than t = kL/c . Electrons,therefore, of progressively higher energy will be confined closer to the inlet, resulting to anenergy-dependent size pipe. This, together with the orientation of the observer, determinesthe time it takes for the emission at a given energy to reach the steady-state.For an observation angle θ = π/
2, there will be no position-dependent delays, given thatthe light path from all parts of the pipe to the observer are equal. Note that if the sourceis moving relativistically with bulk Lorentz factor Γ, θ = π/ θ = 1 / Γ in the 13 –observers’ frame. At ǫ s < ǫ b we expect a practically achromatic increase of the luminosity,reaching a steady state at ∼ kL/c because the electrons responsible for this emission have γ < γ b . At higher energies ǫ s > ǫ b the emission comes from electrons with energy γ > γ b , andthe time to reach steady-state is t = (3 m e c ) / (4 σ τ U B γ ) = (3 m e cB / ⋆ ) / (4 σ τ U B ǫ / s ), where wehave used ǫ s = B ⋆ γ .To verify that the model variability agrees with our analytical predictions, we selectfour synchrotron energies, marked through the four vertical lines in Figure 2. The lowestenergy (dotted vertical line) comes entirely from the ν / tail of the synchrotron emissivity,and it is heavily dominated by the lowest energy electrons with γ = γ min << γ b that haveno time to cool before they escape. The second lower energy (broken vertical line in Figure2) at ǫ min < ǫ < ǫ b is predominately due to electrons with γ min < γ < γ b that also escapebefore they cool appreciably. In the upper panel of Figure 3 we plot the model light curvesof these two low synchrotron energies, using the same line styles. As expected, the two lightcurves are almost indistinguishable, both reaching a steady state at t ∼ L/c (marked by asolid vertical line in the upper panel of Figure 3).The light curves of two more energies, this time with ǫ s > ǫ b , are marked by the dot-dash and triple dot-dash lines in Figure 2 and are plotted in the upper panel of Figure 3using the same line styles. The vertical lines with the same line styles in Figure 3 indicatethe time at which the corresponding light curves are expected to reach a steady state. Ascan be seen, at these times the light curves are at ∼
80% of their steady-state level. This ismainly because electrons continue to radiate at a given energy ǫ s even when their Lorentzfactor drops below ( ǫ s /B ⋆ ) / , as the exponential decay of the emissivity indicates (equation3; e.g. the synchrotron emissivity of an electron with Lorentz factor γ = ( ǫ s /B ⋆ ) / at time t drops ∝ exp[ − ( t/t ) ] from its peak emissivity, requiring t = 2 t to drop by 98%). To veriflythis, we plot in the middle panel of Figure 3 the same light curves, using the δ -functionapproximation for the synchrotron emissivity. As can be seen, the two high energy lightcurves approach the steady-state level significantly closer to the analytically expected time.To evaluate if the light travel effects are properly taken into account, we rotate the pipein such a way that the inlet is closer to the observer ( θ = 0), as is the case for a propagatingshock, that is observed ‘jet on’. In this case, we anticipate the low energy ( ǫ < ǫ b ) variabilitythat requires to fill the entire pipe with plasma, will reach a steady state after an additionaltime L/c , because this is the additional length that variations from the end of the pipe haveto travel to reach the observer. The time it takes, therefore, for the lower energy emissionsteady-state to be reached is ( k + 1) L/c . For higher energy variability ( ǫ > ǫ b ) that reaches asteady-state before time t < kL/c , the additional light travel time required is ( tc/k ) /c = t/k ,where tc/k is the maximum distance from the inlet of the pipe that the contributes to the 14 –energy under consideration. The time it takes, therefore, for the steady-state to be reachedis t (1 + 1 /k ). We demonstrate these considerations in the lower panel of Figure 3. Ascan be seen the time it takes for the low energies ( ǫ < ǫ b ) to reach steady state is now( k + 1) L/c = 3
L/c , and the time it takes for the high energy light curves to reach a steadystate increases by a factor ∼ / t ≈ . L/c for θ = π/ t ≈ . L/c for θ = 0).An analytical result regarding the relation of the synchrotron to the SSC emission, thatwe can test our model against, addresses the relative amplitude of synchrotron and SSCemission. For states that are synchrotron dominated, an increase of the electron injectionnormalization results to a linear increase of the synchrotron emissivity because the syn-chrotron emissivity is proportional to the number of available electrons, and to a quadraticincrease of the SSC luminosity, because the SSC luminosity is proportional to the productof number of electrons times the synchrotron photon energy density, which scales with thenumber of electrons (e.g. Ghisellini, Maraschi, & Dondi 1996). This result holds as long asthe increased injection lasts long enough to occupy the entire volume of the source (in ourcase t var > kL/c ), bringing the source to a new steady state.To verify that our code reproduces this behavior we start from the configuration ofFigure 2, which we observe from an angle θ = π/ t = 3 L/c , longerthan 2
L/c , the time to reach steady state. The middle panel of Figure 4 shows the evolutionof the SED as the flare grows, with the times denoting time since the increased injectionstarted, while the bottom panel shows the evolution of the flare as the flare dies out, startingfrom the time the additional injection is switched off. The characteristic quicker responseof the high energy electrons is apparent. The two solid curves represent the low and highsteady-states. At the upper panel we plot the light curves of the four frequencies markedwith vertical lines at the two bottom panels. They are selected to roughly correspond tooptical, X-ray, GeV and TeV energies, assuming that beaming will increase their observedvalues by a Doppler factor δ ∼ −
40. Note that the synchrotron emission doubles when itreaches its steady-state, while the SSC quadruples, in agreement with the analytical result.
We present here three variability case studies, for which we use as a starting point theconfiguration described in Figure 2, increasing the electron luminosity to L inj = 5 × erg s − to produce a steady-state SED (solid line in the lower panel of Figure 5) that for a 15 –beaming δ ∼ −
40 resembles those produced by flaring TeV blazars. We first study thecase of an increased electron injection that lasts a short fraction of the light crossing time.After the system reaches its steady state, we increase the injected electron luminosity L inj bya factor of 2 for 0 . L/c . At the bottom panel of Figure 5, we plot the SED evolution, whileat the upper panel we plot the light curves that correspond to the four energies denoted byvertical lines at the bottom panel. We note that the flare peaks with no time-delays at allenergies, and then decays with the higher energies of each component decaying first. We alsonote that both the synchrotron and SSC flare have a higher amplitude at higher energies (andtherefore the spectrum hardens as the flux increases). Also, because additional electrons areinjected at all energies, the entire SED responds to the increased injection. The maximumfractional increase of the emission increases with frequency for both the synchrotron and theSSC components. This is because the higher the energy of the electrons required to producea given synchrotron or SSC emission, the shorter their lifetime in the source; an additionalinjection, therefore, for a fixed time (0 . L/c in our case) will increase by a higher factor thenumber of higher energy electrons.To show how the variability event would be seen if we had the ability to resolve thepipe, we plot in the four panels of Figure 6 the luminosity profile along the pipe for the fourenergies we study. In all cases we normalize the luminosity to the steady-state luminosityof the first zone. The lower curve in all cases depicts the steady-state luminosity profile. Asexpected, the high energy synchrotron and SSC emission are confined close to the inlet, whilethe low energy emission of both components extends throughout the source. The decline ofthe steady-state low energy SSC emission along the pipe is due to the fact that this emissionis a convolution of a range of electron energies, and the higher energy electrons are graduallybecoming unavailable away from the inlet. On top of the steady-states, we plot the snapshotluminosity profiles at the times depicted in the lower panel of Figure 5. It can be seen how thepulse is propagating away from the inlet, gradually disappearing due to cooling at the highsynchrotron and SSC energies. At the low synchrotron and SSC energies it can clearly beseen how the pulse is spreading out due to the escape from zone to zone (the n ( γ, t ) /t esc termof the kinetic equation). Note also that while the amplitude of the low energy synchrotronvariation is substantial (starting with an increase by a factor of 2 at the inlet), because thepulse lasts for only 0 . L/c , while low energy synchrotron emission is produced by electronsaccumulating for 2
L/c , the increase of the total low energy synchrotron emission is small(about 5%) as expected, and as can be seen in the upper panel of Figure 5.Another possible variation in the injected electron distribution is an increase in themaximum electron Lorentz factor of the EED, something that can result from a temporaryincrease in the electron acceleration rate. In this case the normalization of the EED remainsconstant and the increase in the injected luminosity depends on the electron index and the 16 –new value of γ max . Using the same configuration as above, we increase the value of γ max by a factor of 5 for the same short time 0 . L/c (for an electron index of p = 1 . γ min = 10 this corresponds to an increase of L inj by ∼ γ ∼ Γ m p /m e , where Γ is the typical bulk Lorentz factor of theflow, to provide the electrons that can be picked up by the Fermi acceleration mechanismfor acceleration to much higher Lorentz factors (e.g. Sikora et al. 2002). Variations in thispre-acceleration mechanism that are not propagated to Fermi acceleration may account forthe variations in the low energy tail of the electron distribution. We simulate this scenarioby injecting an additional low energy EED with the same luminosity as the steady injection,but with γ max = 10 .As can be seen in Figure 8, as soon as the injection starts, the additionally injectedlow energy electrons produce low energy synchrotron emission (dotted line in Figure 8) thatremains at a plateau during the time it takes for the extra injection to transverse the lengthof the pipe, because these low energy electrons do not cool strongly. In the beginning, theselow energy synchrotron photons are produced close to the TeV energy electrons of the steadyflow (due to radiative losses these high energy electrons are found only close to the injectioninlet) and are upscattered by them to TeV energies, producing additional TeV emission whichis manifested as the rising part of the TeV flare (triple dot-dash line). As the event evolves,the additional population of low energy electrons propagates downstream and, due to theincreasing distance between the additionally produced synchrotron photons and the pipe inletwhere the TeV electrons are found, the TeV flare quickly subsides. The behavior of the X-rayemission (broken line) is very interesting: the extra low energy synchrotron photons producedincrease the photon density experienced by the synchrotron X-ray emitting electrons andthis causes additional cooling, slightly reducing the X-ray synchrotron emission. As theproduction of these additional low energy synchrotron photons is displaced downstream,their effect at the neighborhood of the inlet, where the synchrotron emitting electrons arefound, decreases and the X-ray flux returns to its steady state. The behavior of the lowerenergy Gamma-ray emission does not show the sharp decline of the high energy Gammarays. Instead, their light curve shows a gradual decline. This flare is mostly produced 17 –by the additional electrons upscattering a broad energy range of seed photons that has agradually decreasing level with distance from the inlet energy. The model behavior of a TeVflare not accompanied by an X-ray flare is reminiscent of the so-called orphan TeV flares(Krawczynski et al 2004; B la˙zejowski et al. 2005). These are rare flaring states of blazarscharacterized by an increase in the TeV luminosity that is not accompanied by a similarincrease in X-ray energies.
5. Discussion and future work
We presented a multizone code that for the first time takes into account the non-local,time-retarded nature of SSC losses. This code is currently the only multizone model thatincorporates the non-local, time-delayed SSC losses, and as such is uniquely suitable formodeling the results of multiwavelength campaigns at radio, optical, X-ray and γ -ray en-ergies, with the additional constraints in the critical and unexplored for TeV blazars GeVGLAST regime.As we argued, the results of one zone codes for the critical high energy regime of boththe synchrotron and SSC components are problematic, and should not be used to infer thephysical conditions in the source through variability modeling. We described our multizonecode, tested it successfully against known analytical results, and presented a small number ofvariability case studies. The case studies we presented, although based on the same underly-ing steady-state configuration, exhibited very different variability patterns. This means thatdetailed modeling of broadband SEDs and simultaneous multiwavelength variability can beused to infer what is actually the cause of a given observed variability pattern, providingreliable constraints on the particle acceleration taking place. Orphan flares can be repro-duced assuming an increase of the injection of the low energy electrons, but not assumingthe injection of a very high energy electron population, as we also showed analytically.The fact that this plausible variation cannot produce orphan flares significantly narrowsthe parameter space for events that can produce such events, possibly in agreement withtheir observed scarcity.The code we described can run with a typical workstation in a reasonable time of atmost a few minutes at a resolution of ∼
10 bins per decade of observing frequency, ∼ ∼
50 zones. To achieve this we employed a pipegeometry, and adopted an energy conserving δ -function approximation for the SSC emissivity,as well as a step function approach to take into account the change from the Thomson toKlein-Nishina IC scattering cross-sections. Adopting these approximations is problematic for 18 –situations where IC scattering of narrow photon distributions (e.g. line emission from thebroad line region or even a blackbody spectrum characterized by a typical photon energy ǫ )is important. In this case the adoption of the step function cross section description wouldcreate a strong artificial feature on the EED localized at the transition from the Thomson tothe Klein-Nishina regimes at γ ∝ /ǫ , which would then propagate to the emitted spectrathrough the δ -function IC emissivity. For SSC systems, however, where the seed photons arespread over many decades in energy, the resulting spectra are good approximations of thoseproduced using the full expressions for the synchrotron and SSC emissivities as well as thefull Klein-Nishina cross section.Including the above considerations, as well as the processes of synchrotron opacity andpair production through γ -ray absorption within the source, would increase the executiontime up to levels marginally comfortable for the typical workstation. Most probably, suchan extension of the code would require parallelization. A more desirable upgrade of thecode would drop the assumption of no lateral gradients in the plasma characteristics byswitching to a two-dimensional geometry, in which the electron distribution and the SSCphoton energy density are allowed to change laterally to the flow direction. Such consid-erations may be relevant to the recently observed ∼ .
75 days delay between the IR andthe X-ray variability in 3C 273 (McHardy et al. 2007). These authors argued that thedelay may be attributed to the time it takes for the SSC photon energy density to builtup as the SSC photons are transversing the cross section of the flow. This upgrade willscale the computation time roughly by N , increasing it from ∼ few minutes to ∼ sev-eral hours. We note here that our formalism can be extended to treat velocity profiles interm of the decelerating flow (Georganopoulos & Kazanas 2003) or the spine sheath model(Ghisellini, Tavecchio & Chiaberge 2005) that have developed to address the lack of super-luminal motions in TeV blazars (e.g. Piner & Edwards 2004).Another upgrade that can be incorporated in the existing code, this time with a minimalcomputational overhead, is that of a zone for particle acceleration, following the formalismof Kirk, Rieger, & Mastichiadis (1998). In this case, in the first zone of the model, lowenergy electrons will be injected and allowed to accelerate while suffering radiative lossesdue to synchrotron and non-local SSC. These particles will subsequently escape into thepipe and flow downstream. This configuration will require a different numerical schemefor the acceleration zone, since there most particles are advected upward in energy space,but there is a possibility, in a time-dependent scenario, of the highest energy particles beingadvected downward, while the rest of the electrons are still advected upwards. The benefit ofincluding particle acceleration in the code is that it will allow us to study cases of hard lags/counterclockwise loops in the X-ray hardness - X-ray flux diagrams thought to result whenacceleration and loss timescales are comparable. Such a code could be used to model the 19 –observed curved X-ray spectra of high peak frequency blazars in the framework of episodicparticle acceleration (Perlman et al. 2005).We are grateful to the referee for a critical review of the first version of the manuscriptthat made us aware of some important limitations of the first version of our code. Partof this work was done in the context of the senior thesis of Philip Graff at UMBC. Theauthors acknowledge support from NASA LTSA grants REFERENCES
Aharonian, F. A. 2000, New Astronomy, 5, 377Aharonian, F. A., et al. 2005, A&A, 442, 895Aharonian, F. A., et al. 2007, ApJ, 664, L71Begelman, M. C., Fabian, A. C., & Rees, M. J. MNRAS, 384, L19Blandford, R. D. 1978, in Pittsburgh Conference on BL Lac Objects, ed. A. N. Wolfe (Pitts-burgh: Univ. Pittsburgh Press), 328B la˙zejowski, M., Sikora, M., Moderski, R.,& Madejski, G. M. 2000, ApJ, 545, 107B la˙zejowski, M. et. al. 2005, ApJ, 630, 130B¨ottcher, M. 2007, Astrophysics and Space Science, 309, 95Crusius, A. & Schlickeiser, R. 1986, A&A, 164, L16Chang, J. S. & Cooper, G. 1970, J. Comp. Phys., 6, 1Chiaberge, M. & Ghisellini, G. 1999, MNRAS 306, 551D’Arcangelo, F. D. et al. 2007, ApJ, 659, L107Fossati, G. et al. 2000, ApJ, 541, 153Fossati, G. et al. 2008, ApJ, 677, 906Georganopoulos, M. & Marscher, A. P. 1998, ApJ, 506, 621 20 –Georganopoulos, M. & Marscher, A. P. 1998, ApJ, 506, L11Georganopoulos, M. & Kazanas, D. 2003, ApJ, 594, L27Ghisellini, G., Maraschi, L., Dondi, L. 1996, A&AS, 120, 503Ghisellini, G., Tavecchio, F. & Chiaberge, M. 2005, A&A, 432, 401G´omez, J. L., Alberdi, A., Marcaide, J. M., Marscher, A. P. & Travis, J. P. 1994, A&A, 292,33Graff, P. B., Georganopoulos, M., Perlman, E. S., Kazanas, D. 2007, The first GLASTsymposium, AIP Conference Proceedings, 921, 333Harding, A. K. & Lai, D. 2006, Rep. Prog. Phys., 69, 2631Kardashev, N. S. 1962, Soviet Astronomy, 6, 317Kataoka, J., Takahashi, T., Makino, F., Inoue, S., Madejski, G. M., Tashiro, M., Urry, C.M., & Kubo, H. 2000, ApJ, 528, 243Katarzy´nski, K., Ghisellini, G., Tavecchio, F., Gracia, J., & Maraschi, L. 2006, MNRAS,368, L52Kirk, J. G., Rieger, F. M., & Mastichiadis, A. 1998, A&A, 333, 452Krawczynski, H., Coppi, P. S., & Aharonian, F. 2002, MNRAS, 336, 721Krawczynski, H. et. al. 2004, ApJ, 601, 151Maraschi, L., Ghisellini, G., & Celotti, A. 1992, ApJ, 397, L5Maraschi, L. et al. 1999, ApJ, 526, L81Marscher, A. P. et al. 2008, Nature, 452, 966Mastichiadis, A. & Kirk, J. G. 1997, A&A, 320, 19McHardy, I., Lawson, A., Newsam, A. Marscher, A. P., Sokolov, A. S., Urry, C. M. Wehrle,A. E. 2007, MNRAS, 375, 1521Melrose, D. Plasma AStrophysics Vol I, Gordon & Breach, New YorkPerlman, E. S. et al. 2005, ApJ, 625, 727Piner, B. G. & Edwards, P. G. 2004, ApJ, 600, 115 21 –Piner, B. G., Pant, N., & Edwards, P. G. 2008, ApJ, 678, 64Ravasio, M., Tagliaferri, G., Ghisellini, G., & Tavecchio, F. 2004, A&A, 424, 841Sambruna, R. M. et al. 2000, ApJ, 538, 127Sikora, M., Begelman, M. C., & Rees, M. J. 1994, ApJ, 421, 153Sikora, M., B la˙zejowski, M., Moderski, R., Madejski, G. M. 2002, ApJ, 577, 78Sokolov, A., Marscher, A. P., & McHardy, I. M. 2004, ApJ, 613, 725Sokolov, A. & Marscher, A. P. 2005, ApJ, 629, 52Takahashi, T. et al. 1996, ApJ, 470, L89Zhang, Y. H. 2002, MNRAS, 337, 609
This preprint was prepared with the AAS L A TEX macros v5.2.
22 –Fig. 1.— A short variation by a factor of 5 for t inj = t lc / , t lc = R/c , is introduced in theone zone model. The dotted line tracks the synchrotron emission of the low energy electronswith cooling time t cool > t esc , the broken line of electrons with t lc < t cool < t esc , the dash-dot line of electrons with t inj < t cool < t lc , and the solid line of electrons with t cool < t inj .In the upper frame, no light crossing delays are taken into account and this results to theunphysical result of variability times shorter than t lc , for radiation produced by electronswith t cool < t lc . The lower frame depicts the same event with light crossing delays taken intoaccount. In this case, the light crossing time is the smallest observable variability scale. Thetime-integrated emission in the disturbance is the same in both cases, as can be easily seenfor the highest frequency that reaches a plateau in both panels: (5 − × . . − ×
1. 23 –Fig. 2.— Comparison of the steady-steady numerical result with the analytical solution ofa synchrotron dominated configuration. The following parameters have been used: A pipeof length L = 10 cm and radius R = 5 × cm (aspect ratio 10 : 1). The magnetic fieldin the pipe is B = 0 . γ min = 10 , γ max = 2 × , p = 1 .
8. The pipe is split into 40 cylindrical slices of equal height l = L/
40 =2 . × cm, and the escape time is set to k = 2 times the light crossing time. With theseparameters fixed, the ratio of the SSC to synchrotron luminosity increases with increasingelectron injected power, and for L inj = 5 × erg s − the system is synchrotron dominated.At the lower panel we plot the analytical (broken line) and numerical (solid line) steady-stateEED. At the upper panel we plot the analytical synchrotron (broken line) and numerical(solid line) synchrotron and SSC SED. The four vertical lines mark four frequencies for whichwe study their variability in Figure 3 using the same line styles. 24 –Fig. 3.— The light curves of the four frequencies shown by the vertical lines in Figure 2, usingthe same line styles, for two different pipe orientations. The solid vertical lines correspondto the time kL/c required for the low energy ( ǫ < ǫ b ) light curves to reach steady state,while the dot-dash and triple dot-dash vertical lines correspond to the analytical estimatefor the time it takes for the high energy ( ǫ > ǫ b ) light curves to reach steady state. The pipeorientation is θ = π/ θ = 0 for the bottom panel.The middle panel uses a δ -function approximation for the synchrotron emissivity. 25 –Fig. 4.— A quadratic variation: a doubling of the injected electron power for t = 3 L/c ,a time longer than 2
L/c , the time it takes electrons to transverse the pipe. The initialconfiguration is the same as that of Figure 2. The middle panel shows snapshots of the SEDfor a range of times elapsed from the beginning of the additional injection and the bottompanel shows snapshots of the SED evolution after the additonal injection has been switchedoff. The upper panel shows the light curves for the four energies depicted by vertical linesat the lower two panels. 26 –Fig. 5.— Here we use the same steady-state configuration as in Figure 2, except for higherinjected EED luminosity, L inj = 5 × erg s − . After the system reaches the steady-state,we increase the injected luminosity by a factor of 2 for 0 . L/c . The bottom panel shows theevolution of the SED, while the upper panel the light curves of the flare for the four energiesmarked by the vertical lines at the bottom panel. 27 –Fig. 6.— The luminosity profile along the pipe for the four energies whose light curves areplotted in Figure 5, for the same pipe configuration. In each case the lower curve representsthe luminosity profile, normalized to the steady-state luminosity of the first zone. Theprofiles at times t = 0 . , . , . , . , . γ max is increased by a factorof 5 for 0 . L/c . 29 –Fig. 8.— Starting from the same steady state as in Figure 2, we inject for 0 . L/c an additionallow energy EED with the same luminosity as the steady injection, but with γ max = 104