The radio morphology of a spatially resolved SSC model

One of the main, unresolved questions about the nature of quasars is the position of the acceleration site responsible for the highest energies. The attempt to investigate this question in the energy regime with the highest resolution, the radio band, has the downside that no theoretical model exists that can connect these two regimes. The model in this work tries to shrink this gap by extending the general synchrotron self Compton (SSC) model up to length scales in the order of the resolution of radio observations. The resulting spectral energy distributions (SED) show a qualitative improvement in the representation of the radio spectrum. Furthermore the obtained emission morphology shows similar properties to the radio structures observed in jets of quasars. A complete and quantitative connection will however need either much higher numerical effort or an improved methodology.


Introduction
In recent years a diverse set of methods for probing as well as modeling of jets in extragalactic sources has been developed. These methods usually focus on a particular region of either space or energy. Furthermore the different timescales on which these sources change need to be analyzed.
The most detailed view of jets, from kpc to sub-pc scales, can be obtained from radio interferometry with very large baselines (VLBI) [1]. However, the most inner region between the black hole and the radio-core, where the very high energy (VHE) emission is likely to originate from, is not accessible due to self absorption. In the energies above the radio on the other hand, the central engine together with the jet appears as a point source. In contrast to this fact, modeling of the complete spectral energy distribution (SED) above the radio, is very successful by assuming a region of the size in the order of 10 −3 pc. This even holds for very dense energy sampling resulting from multi wavelength observations [2,3]. Especially high peaked BL-Lac sources (HBLs), which show indications of strong Doppler boosting, also show very short variability timescales, supporting the assumption of a very small acceleration and emission region [4,5].
The usually employed synchrotron self Compton models (SSC) [6] are homogeneous in space. Therefore they can not produce radio spectra for the regions and length scales that are observed. On the other hand, large scale jet models depend on (general relativistic) magnetohydrodynamic (MHD) methods. On the one hand, the implementation of acceleration processes is not possible in these a e-mail: srichter@astro.uni-wuerzburg.de b e-mail: felix@fspanier.de models. Hence only synthetic spectra can be used to produce the radiative outputs [7,8]. On the other hand they are scale invariant, so a motivation for a small region that dominates the VHE emission, in order to explain the observed variability, is challenging.
A connection between these distinct regimes could be obtained in the time domain in the form of correlations between variability patterns of the radio fluxes and morphologies and fluxes in bands of higher energies [9][10][11]. However, a theoretical description of this connection does not exist so far.
In section 2 we present a spatially resolved SSC model. It is able to produce non-thermal particle distributions in a small region (i.e. close to a shock) and track it far downstream up to parsec scales. We show in section 3, that this enables both, the interpretation of radio spectra in terms of SSC models, as well as the connection between the region of VHE emission and the radio blobs of VLBI observations. Section 4 will summarize our results. In addition, numerical and conceptional limits of the current model as well as future work will be discussed.

Model
The SSC paradigm is now well established to be responsible for the double humped SED of at least HBLs. However, the shape of the electron distribution is not determined by this paradigm. The power-law form of the SED and the cosmic-ray spectrum suggests to use a broken power-law for the source intrinsic particle spectrum, too. A time dependent treatment, necessary to constrain the model against observed light-curves, should then incorporate a physical sound acceleration mechanism in which the resulting acceleration timescales depend on physical pa-rameters (see e.g. [12]). Moreover, in a spatially resolved model the morphology of the acceleration region has to be well defined.
The most rapid process we know that can produce a power-law distribution is the so called Fermi I acceleration, see e.g. [13]. This process is based on elastic particle scattering in the vicinity of a MHD shock. Acceleration arises then during the isotropization of the particle distribution after the shock crossing. The scattering processes driving this isotropization will take place at various distances from the shock. Hence this process can be accommodated in an inhomogeneous model quite naturally. Any such model including this process should then mimic the Fermi acceleration as detailed as possible to obtain a connection between the variability patterns of the source and the shock parameters.
The here presented model therefore directly connects the shock with the geometry and its discretization. The pitch angle scattering, driving the particle distribution towards isotropy in the rest frame of the ambient plasma, is included stochastically. Since the particles are also advecting through the system and cool simultaneously, this results in a ratio between acceleration and cooling time that depends on the distance to the shock.

Geometry
In this model no test particles are included and only the positions of the gyro centers are computed. These two points limit us to non-relativistic and non-oblique shocks. Consequently, the discretization is carried out along the magnetic field lines, parallel to the shock normal. This is shown for the region around the shock in Fig. 1.  The resulting slices are indexed with i. At each position z i all parameters can in principle be set freely. However, the z dependence of any parameter should be physically motivated. One can then represent the shock by a jump in velocity of the ambient plasma. From the properties of the shock, its velocity V S and compression ratio R, the bulk velocity in each cell is calculated. The downstream velocity behind each shock V P (in units of c), expressed in the upstream frame is From here one can, also for multiple shock scenarios, calculate the bulk velocities for each slice in the frame of the first shock. All other calculations then take place in this frame of reference.
Since the Fermi I process depends on the pitch angle scattering, we can not restrict ourselves to the isotropic approximation. On the other hand, a full discretization of the pitch angle µ = cos(θ) (where θ is the angle between particle momentum and magnetic field) is numerically expensive. Since the pitch angle scattering will produce an isotropic particle distribution on short timescales, it is sufficient to introduce a very crude discretization in the form of two half-spheres: (2) Hence the ratio between n + and n − is a measure for the anisotropy of the distribution. The introduction of a scattering rate between these two quantities (in the rest frame of the bulk plasma), together with the advection between cells, will lead to Fermi-I acceleration. More precisely, pitch angle scattering between uneven distributions will, on average, lead to particle acceleration.
Using the approximation of isotropy per half-space and the particle speed to be V part = 1 (since γ 1) one can compute the average advection speed of the particles, only depending on the bulk plasma in the shock frame V P : Changing the integration boundaries to [−1; 0] will yield V − adv .

Kinetic Equations
All other processes that will influence the energy distribution of the particles are included in the kinetic equation. It is solved time dependent in each cell i and can be derived from the Fokker-Planck equation. Therefor an integration over µ has to be performed. Separating the occurring advection terms yields: The time evolution of all n ± i is computed via the above equation. The processes included in Eq. 4 are the Fermi-II acceleration, synchrotron losses and inverse Compton (IC) scattering. The momentum diffusion coefficient D = v 2 A 9κ depends on the Alfvén speed v A and the parallel diffusion coefficient κ [14]. The synchrotron losses are calculated via where σ T is the Thomsosn cross section [15]. The losses due to the inverse Compton process are calculated with the full Klein-Nishina cross section [16]: The source function S represents the injection at the boundaries of the simulation box. For simplicity, only delta like (in energy) injections at the upstream edge are used. The inverse Compton term (Eq. 6) includes the photon density. Due to this back reaction of the produced photon field a non linear system is created and it becomes necessary to calculate the photon density simultaneously. The time dependent equation for N(ν, t) is The dominant gain ν,sync due to synchrotron emission is computed in the Melrose-approximation [17]: In the same way the synchrotron self absorption coefficient κ ν,S S A is obtained. Similar to Eq. 6 the photon gains and losses due to the IC process are calculated: The catastrophic loss is parameterized by the escape timescale t esc .
The total SED is then calculated similar to the model of Blandford and Königl [18], taking into account light travel times and time dilation.

Results
The spatial resolution of our model offers the possibility to study the emission from different distances from the the shock. The shape of the SED, especially in the radio regime, is closely related to this, as we will show in section 3.3. Before, in section 3.1, the influence of different parameters on the shape of the SED is studied, while in section 3.2 we focus on the adiabatic expansion of the flow behind the shock.

Radio spectral index
The data basis we are starting from is the multi-frequencycampaign of Mkn501 in 2009 [2]. This ensures that the qualitativ study will be performed in a parameter space that is physically relevant. Furthermore a comparison to homogeneous models is possibel. In Fig. 2 two fits are shown. The fitting was performed "by eye". The parameters obtained from these fits are summarized in table 1. They show good agreement with the values found in e.g. [2]. Table 1. Parameters obtained from the fits in Fig. 2 produced with our spatially resolved model. In SSC fits, the radio spectrum is usually ignored. This is due to the fact, that the region observed in the radio is much larger than the one assumed for the higher energies. However, often the radio data is included in such fits, by using a very high injection energy or, in stationary models, a large γ min . The difference between the spectral indeces of the radio in sim1 and sim2 is quite obvious. On the one hand side, this raises the question of the origin of these pre-accelerated particles. On the other hand side, the assumption of a minimum Lorentz factor in a region larger than 10 15 cm would require constant re-acceleration. Otherwise the energies in the particle distribution between γ = 1 and γ min would be populated by synchrotron cooling, when moving away from the shock. Without a physical boundary condition for the emission region, the particles will escape from the shock and cooling will eventually dominate the radio emission and steepen the spectrum due to synchrotron-self-absorption.

Adiabatic expansion
As a more natural way of explaining the observed flat radio spectra we propose an extension of the inhomogeneous SSC model up to the VLBI range (10 18 cm), including the adiabatic cooling of the SSC particle distribution. The simplest geometry is a linear increase of the radius R(z) behind the shock. In addition to the dilution of the particle density the magnetic field will decrease and the particles will undergo adiabatic cooling: Here the zero subscript indicates the values at the shockposition and α denotes the opening angle of the flow. Since we assume a parallel shock and a flow along the direction of the magnetic field we choose m = 2. The radio part of the resulting SEDs is shown in Fig. 3 Figure 3. The radio part of the SED for various simulations of different sizes z and opening angles α. Shown are the runs with α = 0 (no adiabatic expansion) and α = 0.5. The plotted powerlaws reflect the self absorbed limit (s = 3.5), the measured spectral index (s = 1) and the minimum obtained by our simulations (s = 2). The total fluxes are scaled for the sake of readability.
between the size of the simulated region and the spectral index is clearly visible in the case of an adiabatically expanding flow. Although the observed data can't be modeled, this is not expected because of the two orders of magnitude that still separate the sizes of the observed region and the simulation. It can also be observed, that in the current setup an extension beyond a size of z = 10 16 cm only showed marginal improvement. In Fig. 4 of the next section the reason for this behavior will become obvious. Due to the huge numerical costs the question whether this is a systematic effect or is strongly parameter dependent (e.g. the magnetic field structure) has be left to future work at this point.

Emission morphology
The distribution of the emission depending on the distance from the shock is, in principle, directly comparable with VLBI observations. However, resolving four orders of magnitude is numerically quite expensive. Nevertheless one can estimate from Fig. 4 a qualitative pattern.  The photon energies produced by the most energetic particles are confined to a region very close to the shock. This is caused by the rapid synchrotron cooling at these energies. Moving to smaller energies will at first delay the cooling decay, until the self absorbed regime is reached. Here, the fluxes close to the shock are negligible. At larger distances however, where smaller densities and magnetic fields are present, smaller and smaller energies are no longer self absorbed. This results in, with increasing wavelength, growing structures at larger distances from the shock.
In the simulations so far most parts of these structures, even for small radio energies, were already captured within z ≈ 10 16 cm. Consequently, the simulation one order of magnitude larger won't show significant differences in Fig. 3. However, the emerging pattern fits qualitatively well to the core shift (with respect to the wavelength) usually observed in AGN.

Discussion
The numerical model presented in this work extends the canonical SSC model to spatial scales much larger than the region that has to be resolved in order to explain the observed fast variability of the high energy emission. The tracking of the accelerated particle distributions far downstream in different physical environments leads to the conclusion, that a connection between the VHE and radio regime are most likely possible under the assumption of an adiabatically expanding flow. First results presented in this work show an improvement of the fit in the self absorbed regime, when increasing the size of the simulation box up to 10 17 cm. Simultaneously the emerging emission morphology can qualitatively represent radio structures generated from VLBI observations. The quantitative dependence of the emerging structures will be studied in future work. Furthermore, the time correlation between different energy bands usually yield typical time lags between the highest energies and the radio regime. In Fig. 4 there is clearly a potential for time lags between the prompt fluctuations around the shock and the delayed and smeared out variability of the radio structures far downstream. This in turn would allow a prediction for the shock position -and hence the origin of the gamma rays -relative to the maxima of the radio fluxes at different energies. Again, quantitative studies of particular flare scenarios are left to future work.