On the connection of radio and $\gamma$-ray emission in blazars

Blazars are a sub-category of radio-loud active galactic nuclei with relativistic jets pointing towards to the observer. They are well-known for their non-thermal variable emission, which practically extends over the whole electromagnetic spectrum. Despite the plethora of multi-wavelength observations, the issue about the origin of the $\gamma$-ray and radio emission in blazar jets remains unsettled. Here, we construct a parametric leptonic model for studying the connection between the $\gamma$-ray and radio emission in both steady-state and flaring states of blazars. Assuming that relativistic electrons are injected continuously at a fixed distance from the black hole, we numerically study the evolution of their population as it propagates to larger distances while losing energy due to expansion and radiative cooling. In this framework, $\gamma$-ray photons are naturally produced at small distances (e.g. $10^{-3}$ pc) when the electrons are still very energetic, whereas the radio emission is produced at larger distances (e.g. $1$ pc), after the electrons have cooled and the emitting region has become optically thin to synchrotron self-absorption due to expansion. We present preliminary results of our numerical investigation for the steady-state jet emission and the predicted time lags between $\gamma$-rays and radio during flares.


Introduction
Blazars are the most extreme subclass of active galactic nuclei (AGN) having their relativistic jets pointing towards the observer. Their spectral properties are characterized by non-thermal emission over the entire electromagnetic spectrum, rapid variability, high optical polarization and apparent superluminal motion. One of the characteristic features of blazar jet emission is the shape of its spectral energy distribution (SED), which usually has two components: a low-energy component extending from radio to UV/soft X-rays and a high-energy component lying between hard X rays and TeV γ rays.
Multi-wavelength monitoring of blazars offers a unique way of probing the physical conditions in their jets and unveiling the physical processes responsible for their non-thermal variable emission. In recent years, there have been many systematic monitoring programs operating at different energy bands. Their main goal is to investigate the anatomy of blazars by providing high-quality data, probing the rapid variability, and searching for significant correlations between different energy bands, e.g. time lags between radio and γ-ray emission [1][2][3]. Fermi-LAT observations performed over the past 10 years provide crucial spectral and temporal information for hundreds of blazars at energies > 100 MeV [4]. Being an all-sky γ-ray monitor (covering about 20% of the sky [5]), Fermi can detect and report on the flaring activity of blazars, thus often triggering follow-up multi-wavelength observing campaigns. Furthermore, HAWC (High Altitude Water Cherenkov) observatory offers the opportunity to study very high-energy (VHE, E ≥ 100 GeV) flares in survey mode, as it scans 2/3 of the entire sky every day [6]. In an attempt to understand the connection of γ-ray and optical variability in blazars, Robopol (Robotic polarimeter), an optical monitoring program, focuses on the connection between changes in the optical polarization and γ-ray flaring activity in blazars [7].
Multi-year monitoring of blazars at optical and radio frequencies (e.g. BU project [8], F-GAMMA [9]) can provide information about the structure and the kinematics of the jets at sub-parsec to tens of parsec scales. Specifically, at radio many programs play a key role in probing the jet morphology at parcec scales with high resolution radio telescopes and investigating γ-rays and radio variability correlations (e.g. MOJAVE [10], [11], OVRO [12][13][14]).
Despite this multi-year monitoring effort there is still no consensus about the location of the high-energy (X-ray and γ-ray) emission in blazar jets. In particular, the rapid variability, which is an identifying property of blazars, suggests that the non-thermal emission is typically produced in regions of the jet with size 1 light day [e.g. 15,16]. Although electron synchrotron radiation produced in this region can explain the optical-to-X-ray part of the photon spectrum, in most of the cases it cannot account for the jet emission at low frequencies (e.g., ν ≤ 10 10 Hz); in such compact emitting regions, synchrotron radiation is typically self-absorbed [17]. Yet, blazar jets are detected up to GHz frequencies having a power-law spectrum [18]. In order to solve this discrepancy, it has been proposed that the radio emission is in most cases not produced in the same region as the γ-ray photons. The radio-emitting region should be larger in size and less opaque to synchrotron self-absorption [19][20][21][22][23]. The results of the various monitoring programs raise several questions about the dynamical and radiative properties of the jets and their variable radio emission [24][25][26][27]. The localization and the geometry of the emitting source can be estimated by theoretical models in combination with observational data.
Many studies are based on an idealized model of a conical steady radio jet, while the component which is associated with the variability of the source is related to shock waves that are propagating in the jet [28][29][30][31]. Assuming that electrons are accelerated in such shock fronts and are cooled as they move away from it, one gets that the highest frequency synchrotron component is emitted from a thin layer behind the shock, while the lower frequency component is produced in a larger volume, simply because the cooling time for higher-energy electrons is less than for the lower-energy ones [19]. Furthermore, this scenario predicts time lags between high-frequency and low-frequency light curves and a specific evolution of the radio spectrum with time. Building upon these previous studies, [21] proposed a model of a ballistic jet with uniform structure in order to investigate the relationship between low and high energy emission [see also [31][32][33]. In this model, a relativistic electron population is evolved dynamically along the jet, while taking into account the synchrotron and inverse-Compton losses. Then, the synchrotron opacity is self-consistently computed along the line of sight by taking into account the emission at different distances from the jet base.
The goal of this work is to study the relationship between radio and γ-ray activity in blazars and to localize their emission using a simplified framework that allows a wide search of the parameter space. We first construct a phenomenological model for radiative transfer (RT) in a conical flow for computing the non-flaring emission of blazars. More specifically, we apply a time-dependent numerical code that has been originally developed to solve the RT equation (RTE) in a spherical geometry [34] to compute the emission from a conical expanding outflow [see also 35]. To achieve this, we use an expanding spherical volume to solve the RTE [e.g. 36], but "disconnect" the dependence of all other quantities (e.g., magnetic field and particle number density) on the distance from the black hole. In other words, the functional profiles of the magnetic field and particle injection rate are free and not self-consistently computed. We then present a method for studying flaring episodes by changing one or more parameters for a new blob formed at the base of the jet and following this disturbance as it evolves down the jet. In section 2 we outline the model that we adopt and in section 3 we present our preliminary results. We conclude in section 4 with a summary of our results.

Model description and numerical approach
In the present study we construct a parametric leptonic model in order to explain the steady-state spectral energy distribution (SED) of blazars by taking into account radiative and adiabatic losses in an expanding volume [37,38]. The emitting region is moving with a highly relativistic speed v = βc that corresponds to a Lorentz factor Γ = (1 − β 2 ) −1/2 . Then, the Doppler factor of the emitting region is defined δ = [Γ(1 − β cos θ)] −1 , where θ is the angle between the jet axis and our line of sight. Relativistic electrons are injected into the source, that is assumed to be spherical with an initial radius R 0 as measured in its co-moving frame ( see also end of Section 1). The characteristic timescale of the problem is the initial light-crossing time of the blob t cross = R 0 /c, which translates to an observed variability timescale of t var = (1 + Z)R 0 /δc, where Z the redshift; henceforth, we drop all redshift-dependent factors in our analysis. We assume that the emitting region is formed at a distance z 0 from the central engine, with z 0 R 0 , and that its radius increases linearly with time t (as measured in the co-moving frame of the blob) as R(t) = R 0 + u exp t, where u exp is the expansion velocity (u exp < c). The magnetic field strength is parameterized as B = B 0 (R 0 /R) s , where B 0 and s are free parameters of the model. The electrons are confined in the emitting region (i.e., t esc → ∞) and lose energy due to the adiabatic expansion of the source, synchrotron emission, and inverse Compton scattering.
In order to calculate the temporal evolution of the electron distribution function, one has to solve two integro-differential equations, each describing the losses/sinks and injection of relativistic electrons and photons in the emitting region [39]. The kinetic equation of electrons reads: where we used as an independent variable the co-moving source radius R, which is directly related to co-moving time as described above. All terms appearing in the above equation have been transformed accordingly (i.e., the rates are measured per unit distance and not per unit time). The terms A syn , A ICS , A exp are the loss rates for synchrotron emission, inverse Compton scattering, and adiabatic expansion respectively [34]. We model the electron injection rate as: where χ can be positive or negative, and γ min , γ max are, respectively, the minimum and maximum Lorentz factors of the electron distribution. The luminosity of relativistic electrons injected into the blob of radius R can be derived from eq. (2) as L inj e = m e c 2 u exp γ max γ min Q e (γ, R)γdγ. The frequency below which the synchrotron radiation is absorbed can be derived by the condition α ν ssa R ∼ 1 where α ν ssa is the absorption coefficient (e.g. [17]). The synchrotron self-absorption frequency ν ssa is a function of radius through its dependence on the magnetic field and electron number density. In particular, the density of radio-emitting electrons when the blob has radius R can be directly derived from Eq. 1 as C e ≈ t cross u exp V −1 γ 2 where V is the co-moving volume of the blob and the integration over γ's is performed over the relevant range for radio-emitting electrons (the specific values of the integration limits are not relevant for this discussion). We implicitly assumed that the energy loss terms (synchrotron and inverse Compton) for these low-energy electrons are negligible. Using the above expression, one can then calculate the . For all the values of p, s > 0, and χ > 0 we explored, we find that ν ssa decreases for increasing R (in agreement with previous studies). Thus, a blob at small distances from the jet base is optically thick to synchrotron radiation, but it becomes optically thin as it expands and moves to larger distances.
To explain the steady-state blazar emission we first assume that blobs with the same initial properties are continuously produced at a distance z 0 from the central engine. At any given time an observer will receive emission from a series of blobs with different radii from a wide range of distances. This is equivalent to a conical flow with a half-opening angle φ = arctan (R 0 /z 0 ). The distance z traveled by a blob since its "birth", as measured in the black hole's rest frame, is related to its radius R as: z = z 0 + βc(R − R 0 )/Γu exp . Using these assumptions, we integrate along the line of sight the SED in order to reproduce the total steady-state spectrum of the source which is observed. Flaring episodes can be also simulated in this framework, if one blob or a sequence of blobs are injected with different initial properties (e.g. higher value of B 0 or increasing magnetic field with radius (s < 0)) from those producing the steady emission.

Results
In Fig. 1 we present an characteristic example for the steady-state emission of a BL Lac source computed numerically by superimposing the emission of 10 4 blobs that are produced continuously at distance z 0 = 10 −3 pc from the central engine. We follow the evolution of each blob until it reaches the final distance z f inal = 6 × 10 −1 pc. The photon spectra emitted by four indicative plasma blobs that are located at different distance z i from the central engine are overplotted with dashed lines. At the initial time t 0 , the radius of each blob is R 0 = 10 15 cm and its initial magnetic field strength is B 0 = 10 G. The bulk Lorentz and Doppler factors are constant and set to Γ = 5 and δ = 10, respectively. Electrons are initially injected with a luminosity L inj e 0 = 10 42 erg s −1 and with a power-law energy distribution (for numerical values, see figure caption). As the blob radius increases with u exp = 0.3 c, the magnetic field strength and electron injection luminosity decrease as ∝ R −1 (i.e., s = χ = 1). The superposition of emission from multiple blobs produces the steady emission of the jet (black thick line). Figure 2 depicts the light curves at different characteristic energy bands (γ rays, optical and radio) in the case of a flaring episode. As can be seen, there is a time lag between the appearance of the γ ray emission, which is produced immediately, and radio, in which there is a delay due to the synchrotron self-absorption. In the co-moving frame, the time lag could be between few hours to few months depending on the choice of the initial parameters. A first effort to study the effects of the initial parameters on the time lags between the ∼ 3 GHz flux and the 0.4 GeV γ-ray flux is shown in Fig. 3, where we varied systematically one parameter at each time (e.g. R 0 , u exp , L inj e 0 ). Version January 28, 2019 submitted to Galaxies

Summary -Discussion
We have presented a method for calculating the steady-state multi-wavelength blazar emission by approximating the jet flow with a superposition of spherical blobs which, after their production at a fixed distance from the black hole, propagate to larger distances while expanding. We follow the evolution of the relativistic electron population of each blob by taking into account radiative energy losses and adiabatic expansion, and compute the resulting photon emission. We then integrate the emission of all blobs that have propagated at different distances, thus creating an equivalent of a conical outflow, to compute the observed spectrum. This method requires a prescription for the electron injection rate Q e , the magnetic field strength B, and the bulk Lorentz factor Γ of the flow as function of the distance from the base of the jet, or equivalently as a function of the blob radius R.
In the simplest case where Γ is constant, B ∝ R −1 , and Q e ∝ R −1 we find that the bulk of steady-state emission is produced close to z 0 , but the SED is severely self-absorbed at ν ≤ 10 12 GHz. Thus, most of the low-frequency steady-state flux emerges at much larger distances, where both the electron number density and the magnetic field are significantly lower due to the expansion of the source (see Fig. 1), in agreement with previous studies (e.g. [38]). Because of the low photon number density of distant blobs, the absorption of the high-energy photons which are produced at small distances by the more distant ones is negligible. Moreover, the small contribution of the distant blobs to the total steady-state flux suggests that the latter is independent of our choice of the final distance z f inal . Besides the study of the steady-state emission, our approach can be used to model multi-wavelength flares by changing the initial conditions of a blob, namely by mimicking a perturbation in the flow properties at an arbitrary distance from z 0 . The perturbation (or "active" blob) then moves down to the jet and results in a different SED. Figures 2 and 3 show the results of the evolution of such perturbation parametrized by a single blob as it expands. We predict zero or positive time lags, i.e. the γ-rays come first and the radio follow. Our results agree with most of the cases that are observed (e.g. [2]). Alternatively, a flare may be produced by re-acceleration of electrons at a large distance from z 0 . Our approach can also be applied to this scenario by changing the properties (e.g., electron injection rate and magnetic field strength) of a blob after it has reached a certain distance. In the present work we chose to explore cases in which the magnetic field strength and electron injection rate of all blobs (those producing the steady and flaring emission) decrease with distance from the central engine, as we expect in a steady-state jet. In a follow-up extended paper, we will investigate cases with increasing Q e and B when studying flaring episodes.
So far, we have considered a unique and constant Lorentz factor for all the blobs. Nevertheless, our approach can be extended to take into account changes in the bulk Lorentz factor with distance to study the effects of bulk flow acceleration (or deceleration) on both the steady-state and flaring emission [e.g. 29,30]. Although we have presented results for BL Lac objects, we can expand our method by including a varying contribution of external photon fields as the blobs propagate away from the central engine.
In conclusion, we have presented a time-dependent numerical model for computing both the steady-state and flaring multi-wavelength emission from blazars. This model has new ingredients in comparison to previous works: it can take into account synchrotron self-Compton losses, which could be important in some flaring cases, it can be used to investigate flaring events by changing many of the parameters of the blob, and it can treat particle acceleration during flares by including the appropriate terms in the electron equation. Thanks to its parametric nature, the proposed framework allows a wide search of the parameter space and a comparison with previous physical jet models [e.g. 29,31], which will be the topic of a future publication.

Aknowledgments
The authors would like to thank the three anonymous reviewers for their constructive criticism and comments. SB: This research is co-financed by Greece and the European Union (European Social Fund-ESF) through the Operational Programme «Human Resources Development, Education and Lifelong Learning» in the context of the project "Strengthening Human Resources Research Potential via Doctorate Research" (MIS-5000432), implemented by the State Scholarships Foundation (IKY). MP acknowledges support from the L. Spitzer Postdoctoral Fellowship and NASA grant 80NSSC18K1745.