© The Authors 2023

1. Introduction

Blazars belong to the radio-loud subclass of active galactic nuclei (AGN; Blandford & Rees 1978; Blandford & Königl 1979; Urry & Padovani 1995; Blandford et al. 2019; Hardcastle & Croston 2020), with the jet directed along the line of sight of the observer (Dondi & Ghisellini 1995). As a result, the jet radiation is enormously intensified due to the relativistic boosting effect, appearing as one of the most dominant sources for the extragalactic γ-ray sky (Paliya et al. 2019; Hovatta & Lindfors 2019; Bhatta 2022). The typical two-hump structure in its spectral energy distribution (SED) is mainly dominated by two broad non-thermal radiation components (Urry 1999). The low-energy hump is attributed to the synchrotron emission from the relativistic electrons that typically extend from radio to UV or X-ray (Böttcher 2007; Boettcher 2010; Meyer et al. 2012). Often, the presence of different particle energization processes may be accountable for the observed broadening of the low-energy hump beyond the X-ray band (Kirk et al. 1998). The source of the high energy hump is believed to be due to the inverse Compton scattering of low-energy photons and extends all the way to the high-energy end of γ-rays (Liodakis et al. 2018). Furthermore, the production of neutrino emission (IceCube Collaboration 2013) from the blazars suggests pion decay, in addition to proton synchrotron radiation, as the main processes responsible for the high-energy emission (Mannheim 1993; Mücke & Protheroe 2001; Petropoulou et al. 2015).

Based on the optical spectra, the unification scheme of radio-loud AGN classified blazars into two broad categories: flat-spectrum radio quasars (FSRQs) and BL Lac objects (Stickel et al. 1991; Stocke et al. 1991; Urry & Padovani 1995). The Compton dominance is larger in FSRQs due to the presence of external photon sources such as the accretion disk, broad- and narrow-line regions (BLR and NLR respectively), torus, and so on. In these sources, external Compton (EC; Begelman & Sikora 1987; Dermer et al. 1992; Dermer & Schlickeiser 1993; Sikora 1994; Kataoka et al. 1999; Madejski et al. 1999; Błażejowski et al. 2000; Ghisellini & Tavecchio 2009) is the dominant mechanism responsible for the second hump. In certain cases, isotropically present CMB photons are also scattered by the relativistic electrons, giving rise to the IC-CMB process (Böttcher et al. 2008; Meyer et al. 2015; Zacharias & Wagner 2016). However, the BL Lacs are low-power sources and exhibit smaller Compton dominance, suggesting an insignificant contribution from EC. In these sources, synchrotron self-Compton (SSC; Marscher & Gear 1985; Bloom & Marscher 1996) plays a vital role as a possible mechanism for the second hump. The Compton dominance also decreases as the peak frequency of the lower energy hump shifts towards a higher range (Finke 2013; Prandini & Ghisellini 2022). Such features are explained based on the strength of radiative cooling suffered by the emitting electrons for different sources (Ghisellini et al. 1998). In certain blazars, the SED shows unusual characteristics. For example, the multi-wavelength observations of AO 0235+164 show a triple hump structure during its flaring state (Ackermann et al. 2012). Furthermore, the flattening of synchrotron spectra at high energy range is also observed that results in change in the slope of X-ray component in the valley of the SED (Böttcher et al. 2003; Sahakyan & Giommi 2022). However, the origin of such changes in SED from its classical behavior is still unresolved in spite of several propositions that have been put forward to explain such phenomena.

Blazar emission also exhibits multi-timescale variability, with the occasional appearance of quasi-periodicity. Several scenarios have been addressed before as possible explanations for these observed variability timescales. For example, propagation of blobs of plasma through helical magnetic fields (Marscher 2008), presence of shocks in the jets (Valtaoja et al. 1992a,b; Türler et al. 2000; Larionov et al. 2013), geometrical effects related to the viewing angle of the observer with respect to the emission zone (Villata et al. 1998; Raiteri et al. 2017), magnetic reconnection (Ghisellini & Tavecchio 2008; Giannios et al. 2009, 2010; Narayan & Piran 2012; Giannios 2013; Shukla & Mannheim 2020), presence of binary black hole systems (Begelman et al. 1980; Sillanpaa et al. 1988; Gupta et al. 2019), and so on. Recent studies of jets have shown that the optical polarization signatures can be highly variable and correlated with the high-activity state (Abdo et al. 2010; Kiehlmann et al. 2016).

Jet instabilities play a major role in governing the observed signatures of blazar emission. For example, the fluctuations observed in the polarization angle of blazar jets (Zhang et al. 2017) and the quasi-periodic nature of the blazar emission may have a kink origin (Dong et al. 2020). Additionally, the long term variability could also be explained through a kink driven helical jet model (Acharya et al. 2021). Studying the temporal variation of flux provides a broad picture of the source, whereas modeling multi-wavelength spectra is required to understand the substantial and extreme physical conditions within the emission region. In addition, it is also crucial to understand different particle energization mechanisms (such as shock acceleration and magnetic reconnection) and their impact on the light curve and broad-band spectra. In the case of highly magnetized environments, it has been demonstrated that the kink instability generated current sheets might act as possible particle acceleration sites (Bodo et al. 2021; Kadowaki 2021). Some studies also proposed that the blazar flares may be powered by internal shock model (Marscher & Gear 1985; Mimica et al. 2004; Böttcher & Dermer 2010; Moderski et al. 2003; Joshi & Böttcher 2011). Hovatta et al. (2008) investigated the long-term radio variability of a sample of AGN flares and their findings appear to be consistent with the shock-in-jet theory. Fichet de Clairfontaine et al. (2022) have also showed that a strong interaction between the standing shock and moving shock may lead to the generation of flares in the light curve. The above-mentioned studies include several simplifications, such as steady state single zone emission modeling and incorporating simplified or limited high energy emission mechanisms. Moderski et al. (2003) presented a code that simulates the light curves and spectra of blazars during flares by considering a single blob as the emitting region. There have also been approaches to study blazar jet emission by performing 3D numerical simulations with fixed spectra (not accounting for all particle acceleration processes; Dong et al. 2020; Acharya et al. 2021; Kadowaki 2021; Fichet de Clairfontaine et al. 2022). To improve upon these simplifications, we extend our previous work, namely, Acharya et al. (2021) by incorporating consistent emission mechanisms with the hybrid framework of the PLUTO code (Vaidya et al. 2018; Mukherjee et al. 2021).

In particular, we have included a time-dependent multi-zone emission model with suitable radiative and particle acceleration mechanisms (particularly due to shocks). Several numerical studies have combined the macroscopic fluid flow with the microscopic Fokker Planck solver to study shock acceleration. For example, Achterberg & Krulls (1992), Marcowith & Kirk (1999), Wolff & Tautz (2015) adopted a stochastic differential solver to connect these scales. Alternatively, hybrid framework where evolution of a non-thermal electron is followed with the dynamics of the fluid is adopted in several works to study effect of particle acceleration (Micono et al. 1999; Tregillis et al. 2001; Mimica & Aloy 2012; Vaidya et al. 2018; Fromm et al. 2019; Winner et al. 2019; Huber et al. 2021).

Our principal aim is to investigate the impact of shock acceleration on the broadband emission characteristics of blazar jets in the presence of various magnetohydrodynamic (MHD) instabilities. The contents of this paper are divided into two parts. Initially, we focus on the numerical implementation of inverse Compton scattering of external photon fields and understand the contribution of different parameters on the particle spectra and the emissivity. Subsequently, we adopt this implementation to study emission signatures associated with 2D relativistic slab-jet toy model and 3D cylindrical plasma column with different magnetization values.

The paper is structured as follows. In Sect. 2, we describe the numerical implementation of the external Compton process for a mono-directional source of seed photons. This section includes the calculation of energy loss rate and external Compton scattering emissivity. Section 3 is devoted to the results obtained from a 2D relativistic slab jet simulation along with the numerical setup and emission modeling approach. In Sect. 4, we explain the multi-wavelength nature of a relativistic plasma column that has undergone kink instability and present a comparative analysis among different parameters responsible for the highly energetic external Compton process. Finally, we discuss all the results and summarize our current findings in Sects. 5 and 6, respectively.

2. Implementation of external Compton mechanism

We adopted the PLUTO code (Mignone et al. 2007) for all studies in the current work; specifically, we used the hybrid framework in the PLUTO code (Vaidya et al. 2018; Mukherjee et al. 2021). This framework is developed to model the non-thermal spectral signatures of the macro-particles in the relativistic MHD flow. The macro-particles are assumed to be an ensemble of leptons with a finite energy distribution. Based on the fluid conditions, the spectral distribution of each particle is updated with time. This enables us to account for several physical processes and emission mechanisms, such as synchrotron and inverse Compton (IC) scattering. The current framework already has a numerical scheme to account for energy losses due to synchrotron and IC-CMB (IC scattering of CMB photons) processes. In this work, we have additionally incorporated the IC scattering process, where the origin of the target photon field considered for the IC scattering is external to the jet, such as the BLR region, accretion disc, and so on. Such an emission mechanism is typically known as the external Compton (EC) process. In particular, we have incorporated the EC process, considering a mono-directional photon field and all the formalisms provided in this regard are taken from Khangulyan et al. (2014).

2.1. Calculation of loss rate

A semi-analytical approach is considered by Vaidya et al. (2018) for the evolution of the spectral distribution of the macro-particles using a Lagrangian scheme. For this purpose, they solve a characteristic equation to calculate the energy loss rate due to the above-mentioned physical processes and emission mechanisms. We extend the energy loss rate equation with an additional term that accounts for the loss due to the EC process:

$$\begin{array}{c}\hfill \frac{\mathrm{d}E}{\mathrm{d}{\tau }_{\mathrm{p}}}=-c_1({\tau }_{\mathrm{p}})E-c_2({\tau }_{\mathrm{p}})E^2-c_{3}Ef(E)\equiv \dot{E},\end{array}$$(1)

where the first term in Eq. (1) represents the loss due to adiabatic expansion, the second term represents the loss due to both synchrotron and IC-CMB (under the Thomson limit), and the third term corresponds to the loss due to the EC process. Also, τ_p is the proper time measured in the jet co-moving frame. The constants associated with all these loss terms are given below:

$$\begin{array}{cc}\hfill c_1& =\frac{{\mathrm{\nabla }}_{\mu }{u}^{\mu }}{3}\hfill \\ \hfill c_2& =\frac{4{\sigma }_{\mathrm{T}}c{\beta }^2}{3m_{\mathrm{e}}^2{c}^4}[U_{\mathrm{B}}+U_{\mathrm{rad}}(E_{\mathrm{ph}})]\hfill \\ \hfill c_3& =\frac{2r_0^2{k}_{\mathrm{B}}^3{T}_{\mathrm{cmv}}^3{\kappa }_{\mathrm{cmv}}}{\pi c^2{ħ}^3}·\hfill \end{array}$$(2)

In the above equations, c₁ can be calculated from the mass conservation equation as $\frac{{\mathrm{\nabla }}_{\mu }{u}^{\mu }}{3}=\frac{1}{3\rho }\frac{\mathrm{d}\rho }{\mathrm{d}{\tau }_{\mathrm{p}}}$ with u^μ being the bulk four-velocity and ρ is the density of the fluid. Here, β and m_e are the velocity (in the units of speed of light, c) and mass of the electrons, respectively, and σ_T is the Thomson cross-section. The quantities U_B and U_rad are magnetic and the radiation energy densities, respectively, and E_ph is the energy of the incident CMB photon, while r₀ is the electron classical radius, κ_cmv and T_cmv are the dilution factor and temperature of the photon field in the comoving frame of the jet. Boltzmann constant and Planck’s constant are denoted in usual manner as k_B and ℏ = h/2π, respectively. It is also important to take notice of the fact that the temperature and the dilution factor should be corrected for the relativistic bulk motion of the emitting region in the following way (Khangulyan et al. 2014):

$$\begin{array}{c}\hfill T_{\mathrm{cmv}}=\frac{T}{D_{\mathrm{PH}}},{\kappa }_{\mathrm{cmv}}=D_{\mathrm{PH}}^2\kappa ,\end{array}$$(3)

where T and κ represent the photon temperature and dilution factor as measured in the lab frame. These quantities are set as input parameters in our model. D_PH = [Γ(1 − β cos χ)]⁻¹ is the Doppler factor between the source of seed photons and the EC emitting (jet) region, where χ is the angle made by the velocity vector of the electron with the photon’s direction. For the present work, we have considered mono-directional photon field such that it lies beneath the EC emitting (jet) region and along the direction of jet motion (such that χ ∼ 0). Inverse Compton scattering is therefore implemented by assuming a tail-on collision where the electron’s velocity vector makes a zero degree angle with respect to the direction of the photon.

The dilution factor for a mono-directional photon field can be approximately represented as the following:

$$\begin{array}{c}\hfill \kappa =\frac{\mathrm{\Delta }\mathrm{\Omega }}{4\pi },\end{array}$$(4)

where ΔΩ is the solid angle of the target photon field as observed from the EC emitting region. In our work, ΔΩ ≪ 1 since the emitting region is considered tens of parsec away from the external photon field and can therefore be further expressed as:

$$\begin{array}{c}\hfill \kappa ={\left(\frac{r_{\mathrm{PH}}}{2r_{\mathrm{EC}}}\right)}^2,\end{array}$$(5)

where r_PH is the radius of the target photon field and r_EC is the distance between the source of the target photon field and the EC emitting region. The cartoon representation of this kind of formalism is shown in Fig. 1. There is a super massive black hole (SMBH) at the centre of the system surrounded by the BLR region, a dusty torus with a relativistic jet lying perpendicular to the plane of the accretion disk. As the emitting region is far from the central zone, it can be considered as a point source when viewed from the emitting region.

Fig. 1. Cartoon representation of emitting region and external photon fields with the seed photon’s direction. The size and distance of different components are not shown to scale.

The energy dependence term in Eq. (1), that is, the energy-loss rate of an electron due to the EC process in the presence of an external photon field in the jet comoving frame is adopted from Khangulyan et al. (2014) and is given by

$$\begin{array}{c}\hfill f(ϵ)=\frac{c_{\mathrm{e}}log(1+0.722\times \frac{ϵ}{c_{\mathrm{e}}})}{1+(c_{\mathrm{e}}\times \frac{ϵ}{0.822})},\end{array}$$(6)

where c_e = 4.62 and ϵ = 4γΘ, $\gamma =\frac{E}{m_{\mathrm{e}}{c}^2}$, and $\mathrm{\Theta }=\frac{k_{\mathrm{B}}{T}_{\mathrm{cmv}}}{m_{\mathrm{e}}{c}^2}$ are the energy of the electron and temperature of the photon field in the jet comoving frame, respectively, given in the units of m_ec². Such approximation for estimation of f(ϵ) is applicable in both Thomson and Klien Nishina limits and provides accuracy of an order of 1%. Besides, as was already noted, the dilution factor takes care of the angular distribution of seed photon field, and the temperature of the black body inherently is dependent on the energy density of the photon field. These presumptions are addressed by the analytical expression that is presented here. In this work, we solve Eq. (1) numerically by implementing the fourth order Runge-Kutta method. The particle spectra and the energy distribution in each bin are updated using the formalism provided in Vaidya et al. (2018). A brief discussion on the interplay among the multiple radiation mechanisms in the cooling process is discussed in Appendix A.1.

2.2. Calculation of the external Compton emissivity

In this section, we describe the quantification of emissivity due to EC scattering due to the interaction of relativistic electrons with a given radiation field. In the present work, we assume that the black body photon field is mono-directional and adopted an expression of the interaction rate with relativistic electrons in the co-moving jet frame.

For the case of mono-directional black body radiation field, the interaction rate in the co-moving frame is given by (Eq. (11) of Khangulyan et al. 2014):

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{ph}}}{\mathrm{d}{\omega }_{\mathrm{sc}}\mathrm{d}t}=\frac{2r_0^2{m}_{\mathrm{e}}^3{c}^4{\kappa }_{\mathrm{cmv}}{\mathrm{\Theta }}^2}{\pi {ħ}^3{\gamma }^2}[\frac{q^2}{2(1-q)}{F}_1(x_0)+F_2(x_0)],\end{array}$$(7)

where ω_sc = qγ is the up-scattered photon energy in the units of m_ec² and x₀ = q/(1 − q)ϵ_ϑ, with ϵ_ϑ = 2γ Θ(1 − cos ϑ_e), $cos{\vartheta }_{\mathrm{e}}=\frac{cos{\theta }_{\mathrm{obs}}-\beta }{1-\beta cos{\theta }_{\mathrm{obs}}}$ is the electron’s direction, and θ_obs is the observing angle. The two analytic functions F₁ and F₂ are estimated with 1% accuracy, as demonstrated in Khangulyan et al. (2014).

We further assume that the relativistic electrons within each macro-particle are isotropically distributed in the co-moving jet frame. The assumption of such isotropic electron distribution is indeed a simple, zeroth-order approximation and it allows us to simply determine the electron direction; it was previously used in Dermer (1995), Georganopoulos et al. (2001). We can therefore estimate the total emissivity of the electron population by convolving the above interaction rate with electron distribution as follows:

$$\begin{array}{c}\hfill {\nu }_{\mathrm{sc}}{j}_{\mathrm{EC}}^{\prime }({\nu }_{\mathrm{sc}})={\displaystyle {\int }_{{\gamma }_{\min }(t)}^{{\gamma }_{\max }(t)}\left[\frac{\mathrm{d}{N}_{\mathrm{ph}}}{\mathrm{d}{\omega }_{\mathrm{sc}}\mathrm{d}t}\right]{\omega }_{\mathrm{sc}}{m}_{\mathrm{e}}{c}^{2}N(\gamma )\mathrm{d}\gamma .}\end{array}$$(8)

In the above equation, the scattered frequency, ${\omega }_{\mathrm{sc}}=\frac{h{\nu }_{\mathrm{sc}}}{m_{\mathrm{e}}{c}^2}$, where ν_sc can be obtained from the observed frequency, ${\nu }_{\mathrm{sc}}=\frac{{\nu }_{\mathrm{obs}}}{D}$, with D = [Γ(1 − β cos θ_obs)]⁻¹ is the Doppler boosting factor. The limits to the integration are the minimum and maximum value of the electron energy in terms of m_ec² and are time dependent as the electron energy is evolving following the radiative loss Eq. (1).

The emissivity in the observer frame in the units of erg cm⁻³ s⁻¹ Hz⁻¹ str⁻¹ is given by

$$\begin{array}{c}\hfill j_{\mathrm{EC}}^{\mathrm{obs}}({\nu }_{\mathrm{obs}})=D^2{j}_{\mathrm{EC}}^{\prime }({\nu }_{\mathrm{sc}}).\end{array}$$(9)

Finally, we can deposit the EC emissivity obtained for each macro-particle on to the grid cells so as to give the grid distribution of $j_{\mathrm{EC}}^{\mathrm{obs}}({\nu }_{\mathrm{obs}},\mathit{r})$ and then can obtain specific intensity maps by integrating along the line of sight.

3. Slab jet simulation

To understand the fundamental impact of high-energy emission mechanisms such as the EC process on the energy spectra and emission features, we first focus on a two-dimensional relativistic slab jet. For that purpose (in Sect. 3.1), we discuss the numerical setup required to simulate a 2D relativistic slab jet along with emission modeling particulars. Section 3.2 discusses the results obtained from our simulations in the context of multi-band emission. Furthermore, a brief discussion on the evolution of particle spectra due to several radiative processes, in addition to the validation of our numerical algorithm, is provided in Appendix A.2.

3.1. Model setup

In this work, the numerical simulations are carried out using the relativistic magnetohydrodynamic (RMHD) module of the PLUTO code (Mignone et al. 2007). The numerical setup is the same as described in Vaidya et al. (2018). The summary of the computational details and initial value of the parameters are given in a tabular form (see Table 1). The simulation is carried out in a Cartesian domain having a dimension of x = (0, L) and y = ( − L/2, L/2). The jet of a width, l, is under-dense in comparison to the ambient medium with a density ratio, η. Furthermore, an initial uniform magnetic field is applied along the x-direction corresponding to a plasma beta of value 10³. The flow has a velocity in the same direction where the ambient medium is static. The initial values of the magnetic field strength and Lorentz factor is given in Table 1. Following Bodo et al. (1995), we apply a perturbation in the y component of the velocity that leads to the development of Kelvin-Helmholtz instability (KHI) and consequently generates shocks at the jet boundary or at the vortices. The total computational time of the simulation is 200 t_sc, which is equivalent to a duration of 0.13 Myr in physical units with t_sc = 6.52 × 10² yr. Furthermore, the unit density on the central axis of the jet and the unit length scale is chosen to be 1.661 × 10⁻²⁸ gm cm⁻³ and L = 2000π pc, respectively.

We performed two different simulations: case 1 – without the inclusion of energy loss due to EC; case 2 – including the energy loss due to EC. These simulations are used to present a comparative analysis among the energy loss terms given in Eq. (1) and the corresponding results are given in Sect. 3.2 and Appendix A.2. To study the effects of different energy-loss mechanisms and their impacts on the emission signatures, we introduced 2 × 10⁵ number of macro-particles at the initial time of the simulation. The distribution of the Lagrangian particles allows a complete sampling of the system, as shown in Fig. 2. Here, the particles are coloured based on their unique identity, ranging from 1 to the maximum number of macro-particles. The background on which these particles are over-plotted represents the fluid density in grayscale. The left and right panels of Fig. 2 are the snapshots taken at times t/t_sc = 150 and 200, respectively. The combined evolution of Lagrangian macro-particles and the fluid density from our 2D runs clearly shows that the slab jet is sufficiently sampled during the evolution. Further, one can observe that the vortices of the KHI are responsible for mixing the particles within the slab jet as it evolves. An initial spectra with a power law index of p = 6 is considered for each macro-particle, with energy cutoffs of γ_min = 10² and γ_max = 10⁸ distributed over 256 bins. It should be noted that the choice of particle spectra is arbitrary. However, since the spectra will eventually flatten as a result of the particles being shocked, we have initially chosen a steeper index. The synchrotron emissivity for an initial power law particle spectra is calculated using Eq. (37) of Vaidya et al. (2018) and the EC emissivity is calculated using Eq. (8) of Sect. 2.2. All the emissions are obtained at the time t/t_sc = 200 (particle distribution is shown in the right panel of Fig. 2) with an observer making a 5° angle with respect to the z-axis (pointing out of the plane).

Fig. 2. 2D representation of particle distribution at t/tsc = 150 (left) and 200 (right) overlaid with the fluid density (grey color). The colour bars show the fluid density and particle unique identities.

3.2. Effect on multi-band emission maps

Another way of studying the effect of energy loss due to the EC process is to estimate the emissivity at different frequencies. The temperature of the seed photon field is assumed to be 1000 K (corresponding to typical temperature of IR photons from the torus Donea & Protheroe 2003; Malmrose et al. 2011; Oyabu et al. 2017) for the calculation of EC emissivity, and the dilution factor value is estimated to be κ = 10⁻⁷ by considering the r_PH ≈ 0.1 pc (representative scales associated with inner structure of torus Urry & Padovani 1995) and r_EC > 10 pc. The top and bottom panels of Fig. 3 show the emission maps for case 1 and case 2 in optical R-band (left) and 10¹⁷ Hz (equivalent to 0.3 keV; right) respectively. The emissions shown at both frequencies are normalized to their individual maximum values and given in the units of ergs s⁻¹ cm⁻³ Hz⁻¹ str⁻¹. The maximum values of emissivities for R-band and 10¹⁷ Hz are 1.18 × 10⁻⁴⁰ and 1.74 × 10⁻⁴⁷ (top) and 3.44 × 10⁻⁴¹ and 6.25 × 10⁻⁴¹ (bottom) respectively. Such a comparative analysis of emission due to different radiative processes exhibits different structures observed in the jet. The emission maps for case 1 and case 2 show significant differences at both frequencies (top and bottom left-hand panels). For case 1, the synchrotron emission is mostly coming from the shocked region. However, for case 2, an extended emission is observed at both R-band and 10¹⁷ Hz, with the inclusion of loss due to EC compared to the case with loss due to only synchrotron and IC-CMB. This suggests that for the chosen set of values of photon field temperature and κ, the addition of loss due to the EC process significantly affects the particle spectra. It can be seen from the time evolution plot of particle spectra for a single test particle (see the bottom panel of Fig. A.2) which in turn reflects on the emission signatures. The lower energetic particles lose energy significantly because of EC loss compared to the high energy particles because of synchrotron loss. These particles up-scatter the lower energetic photons and are responsible for the observed extended emission. In case 1, however, there are no such electrons to emit optical and X-rays throughout the jet, hence, we only observe emission at the shocked region. We note that for the scenarios with a lower photon field temperature or lower κ values compared to the considered ones, the effect of loss due to EC will be reduced. As a result, the EC emission were then observed at 10¹⁷ Hz and at even higher energies.

Fig. 3. Emissivity slices of slab jet for case-1 and case-2 at an observing frequencies ν = R-band and 1017 Hz at t/tsc = 200. Here, the emissivity values are normalized to their individual maximum values.

Additionally, we estimated the emission at multiple frequencies, starting from radio to γ-rays. The top panel of Fig. 4 shows emissivity maps obtained at 1.4 GHz, 43 GHz, and optical R-band; in the bottom panels, we show the emission maps are obtained at ν = 10¹⁷ Hz, 10²⁰ Hz, and 10²³ Hz, respectively. The emissivity values are normalized to their individual maximum values, showing different features at different frequencies and given in the units of ergs s⁻¹ cm⁻³ Hz⁻¹ str⁻¹. The maximum values at 1.4 GHz, 43 GHz, R-band, 10¹⁷ Hz, 10²⁰ Hz, and 10²³ Hz are 1.59 × 10⁻³⁵, 5.5 × 10⁻³⁹, 3.44 × 10⁻⁴¹, 6.25 × 10⁻⁴¹, 1.07 × 10⁻⁵¹, and 1.29 × 10⁻⁶⁰ respectively. As expected, in the radio bands, the synchrotron emission gets diminished with an increase in the observing frequencies, and the features get enhanced at 43 GHz due to the lower normalization value. In the optical and 10¹⁷ Hz (X-ray), an extended emission is observed throughout the jet in addition to the shocked emission at the sheared regions. The emission weakens due to radiative cooling and additional features are visible at higher energies. In the EC scattering, the lower energetic photons gain energy and get up-scattered, where the frequency of the photon after up-scattering is given by ν_up ≈ γ²ν₀ (Longair 2011). In our study, the considered photon field temperature of T = 1000 K corresponds to IR photons with ν₀ ≈ 10¹³ Hz. According to the recipe given in Longair (2011), the IR photons would require electrons of Lorentz factor γ ∼ 10¹ − 10³ to obtain emission in the R-band and at 10¹⁷ Hz. Similarly, electrons with γ ∼ 10⁴ and γ ∼ 10⁵ − 10⁶ are required to obtain emissions at 10²⁰ Hz and 10²³ Hz, respectively. It is expected to have a higher number of lower energetic electrons since initially we considered a power-law particle spectra from γ_min = 10² to γ_max = 10⁸. In addition, due to the strong effect of loss due to EC, many particles would lose energy and reach up to γ_min = 10¹. Such a distribution for the lower energy electrons is responsible for the extended optical and X-ray emission. As the number of particles is reduced with the increase in γ of the electrons, the emission at higher frequencies is also reduced.

Fig. 4. Emissivity slices of slab jet for case 2 for different observing frequencies at t/tsc = 200. Here, the emissivity values are normalized to their individual maximum values.

4. Plasma column simulation

Blazar jets exhibit multi-timescale variability and emit non-thermal emissions covering the entire gamut of the electromagnetic spectrum. These jets are highly magnetized and are prone to undergo several MHD instabilities during their propagation in space and could possibly trigger jet radiation and particle acceleration. In Acharya et al. (2021), we investigated the physical configurations preferable for the formation of kink mode instability by performing relativistic MHD simulations of a plasma column using the PLUTO code. The plasma column depicts a portion of an AGN jet about a few tens away from the central engine and hence magnetically dominated and relativistic in nature. Additionally, connecting the dynamics of the plasma column with its emission features, we found a correlated trend between the growth rate of kink mode instability with its flux variability, obtained from the simulated light curves. We would like to stress that in our previous work, we focused only on the optical band and have estimated the synchrotron emission by using a static particle spectra. However, in this work, we have improved our emission modeling approach by using an evolving particle spectra and including a high energy emission mechanism for a better understanding of multi-wavelength emission properties.

4.1. Model setup

The simulations are initialized in the computational domain of size 8 R_j × 8 R_j × 12 R_j with R_j = 0.5 being the radius of the cylindrical plasma column. A force-free initial magnetic field profile is chosen for the equilibrium condition (Mizuno et al. 2011; Anjiri et al. 2014) and is given by

$$\begin{array}{c}\hfill B_z\frac{\mathrm{d}{B}_z}{\mathrm{d}R}+\frac{B_{\varphi }}{R}\frac{d}{\mathrm{d}R}(R{B}_{\varphi })=0,\end{array}$$(10)

where, $R=\sqrt{x^2+{\mathit{y}}^2}$ is the radial position in the cylindrical coordinate system. Here, the magnetic field in the radial direction is B_r = 0, with B_z and B_ϕ being the poloidal and toroidal magnetic field components. The column has a flow velocity in the z-direction, given by a bulk Lorentz factor, Γ_z = 5, where the ambient medium is static. The initial equilibrium of the system is perturbed by a velocity provided in the radial direction. We note that the perturbation is provided in such a way that the number of kink that fit into the simulation box is n = 4. Following Acharya et al. (2021), the on-axis column magnetization (magnetic-to-matter energy density) in the relativistic form is given by ${\sigma }_0=\frac{B_0^2}{{\rho }_{\mathrm{c}}}$. Here, B₀ and ρ_c = 1.0 are the magnitudes of the magnetic field strength and density on the axis given in the non-dimensional units, respectively. To study the effect of magnetization on the emission features, we chose two values of σ₀. The other detailed description of the initial conditions and numerical methodology is provided in Acharya et al. (2021). Here, the dimensionless quantities can be scaled with appropriate physical units, relevant for the present study. For the relativistic module, the unit velocity is equal to the speed of light c = 2.998 × 10¹⁰ cm s⁻¹, the unit length is chosen to be R_sc = 0.1 pc and the unit density is ρ_sc = 1.673 × 10⁻²⁴ gm cm⁻³. The above choices set the time and the magnetic field to be in units of t_sc = 0.32 yr and B_sc = 1.374 × 10⁻¹ Gauss, respectively.

To understand the effect of the dynamics of the system on the non-thermal emission, we initialized all the runs with 3 × 10⁵ number of Lagrangian macro-particles using the hybrid framework of PLUTO code. This methodology takes into account the effects of microphysical processes on the distribution functions of emitting particles and, subsequently, on emissivities. To model the non-thermal emission, similarly to the slab jet problem, we considered a power-law particle spectra with an initial power law index p = 6 with energy bounds γ_min = 10² and γ_max = 10⁸ distributed over 256 bins with equal bin widths on logarithmic scale. The synchrotron and EC emissivities are obtained using the equations mentioned in Sect. 3.1, while considering an observer making a 5° angle with respect to the axis of the column (jet). In this work, we have studied two reference cases with σ₀ = 10 and 1, named Ref_s10 and Ref_s1. In addition, for comparative analysis, different seed photon field temperatures and κ values are considered, and the details of these simulation setups are given in Table 2.

4.2. Results

In this section, we mainly focused on the results obtained from the Ref_s10 and Ref_s1 cases and discussed their multi-wavelength properties (simulated light curve and SED). The variability and correlation studies provide useful information regarding the location and size of different emitting regions and sources of high energy emission in differently magnetized environments.

4.2.1. Emission maps

The contribution of loss due to EC on the emission features has already been shown in Sect. 3.2 in the context of a 2D relativistic slab jet. Here, in Fig. 5, we have shown the X − Z cuts of normalized emissivity slices at Y = L_y/2 for Ref_s10 case at t/t_sc = 26 (top) and 70 (bottom). At both time stamps, in radio bands (see panels a, b, f, and g), strong emission is coming from the boundary in addition to the extended emission observed throughout the column. With the onset of the instability at t/t_sc = 26, there are formation of shocks. As a result, we observed localized synchrotron emission near the axis of the column in R-band. However, at this early evolution time of the instability, the loss rate is not sufficient to exhibit emission at 10¹⁷ Hz. At higher energies (see panel e), EC is the dominant process and exhibits an extended emission throughout the column. Due to the high growth rate of the instability in Ref_s10, the emissivity maps show prominent structures of the formation of kink at t/t_sc = 70. At this time stamp, in the optical band, the localized synchrotron emission is mostly observed at the shocked positions or at the kinked portion of the column. At 10¹⁷ Hz, an extended EC emission is observed throughout the jet from the contribution of lower energetic particles, as a result of strong EC cooling. Additionally, there is a contribution of localized synchrotron emission from the shocked positions.

Fig. 5. X − Z cuts of emissivity slices for Ref_s10 case at observing frequencies ν = 1.4 GHz (a), 43 GHz (b), R-band (c), 1017 Hz (d) and 1020 Hz (e) at time stamps t/tsc = 26 (top) and 70 (bottom). The emissivity values are normalized to their maximum values and given in the units of erg s−1 cm−3 Hz−1 str−1. The normalized values of 3.42 × 10−22, 4.23 × 10−22, 3.94 × 10−25, 8.31 × 10−27, 3.4 × 10−35, 2.83 × 10−22, 1.74 × 10−22, 4.76 × 10−26, 6.71 × 10−29, and 3.4 × 10−38 for panels a, b, c, d, e, f, g, h, i, and j, respectively.

In an environment of low magnetization (i.e., the Ref_s1 case) due to the trans-Alfvénic nature of the flow, a mixing of kink and KHI is expected. However, as a consequence of lower magnetic field strength, the growth of the instability is not enough to exhibit notable morphological changes (emission maps are not shown here). The synchrotron emission gets diminished while observing at radio frequencies due to the radiative cooling effect. In the optical band and 10¹⁷ Hz, the very few localized shocked particles emit synchrotron emission in addition to the broad or extended EC emission throughout the column. A comparably lower emission is observed at 10²⁰ Hz as a result of EC cooling.

4.2.2. Multi-wavelength light curve

One of the ways to understand the emission characteristics is to estimate the total integrated flux and study its time-varying properties. All relativistic effects, such as relativistic boosting and light travel effects, have been taken into account while estimating the integrated flux. In particular, the emissivities are obtained in the observer frame by multiplying the comoving emissivity by the square of the Doopler boosting factor (e.g., see Eq. (9)). Furthermore, we have used a realistic slow-light approximation (Bronzwaer et al. 2018) to collect information about the light as it travels through the grid cells of the simulation box. We have shown the simulated multi-wavelength light curve for the Ref_s10 case in Fig. 6, where the fluxes are normalized to their maximum values and provided in the units of erg s⁻¹ cm⁻² Hz⁻¹. The light curve is shown from t/t_sc = 20 to 80, corresponding to a period of ≈20 yr in physical units. During this period, the instability is in its evolving state and, hence, it is favourable to study the associated emission signatures during this phase. The yellow highlighted or shaded region indicates a period when the light curve shows distinctive features and mainly corresponds to the linear growth phase of the instability. Later on, during the non-linear growth of the instability, the magnetic field structure becomes chaotic and a decaying phase of the light curve is observed. The black vertical dotted lines are plotted to indicate different activity states of the system and are labelled as state 1, state 2, state 3, and state 4. In state 1, at t/t_sc = 26, an enhancement of the flux can be seen in the given frequencies shown in Fig. 6. In state 2, at t/t_sc = 33, a less active state is observed in all energy bands compared to state 1. In state 3, at t/t_sc = 39, a semi-harmonized behavior along with flux enhancement is observed, suggesting a moderately active state of the system. Lastly, in state 4, at t/t_sc = 70, a comparably low-activity state is observed since not much variation in flux can be noticed at that particular time. The spectral behavior of the system at these particular times is also given in Sect. 4.2.4. In radio bands, a moderately variable emission is seen during the yellow shaded duration. However, in the optical R-band and 10¹⁷ Hz, a characteristic transient feature is seen at t/t_sc = 26, followed by a sharp decay in the light curve. Further, at t/t_sc = 39, few particles encounter shocks again leading to the high flux state in R-band and 10¹⁷ Hz compared to the state-2 (t/t_sc = 33) where the previously shocked or emitting particles were in a cooling phase. Such transient characteristics are due to the contribution of localized synchrotron emission as a result of the generation of shocks. We note that the flux values are very low at the higher energies compared to the lower energies and the certain extreme jump in flux in R-band and 10¹⁷ Hz is due to the few number of particles encountering localized shocks. A decaying nature of the light curve is observed at 10²¹ Hz (γ-ray).

Fig. 6. Simulated light curve for the Ref_s10 case for an observer making 5° angle with the axis of the column. The observing frequencies are ν = 1.4 GHz, 43 GHz, 1.4 × 1014 Hz (R-band), 1017 Hz, and 1020 Hz, respectively, from top to bottom. The vertical black dotted lines and the shaded regions correspond to different activity state of the system (discussed in Sect. 4.2.2). The flux values from top to bottom are normalized to their maximum values 2.64 × 10−17, 2.52 × 10−17, 2.66 × 10−21, 5.13 × 10−24, and 1.54 × 10−31, respectively, and given in units of ergs s−1 cm−2 Hz−1.

Similarly, we have shown the simulated light curve for the lower magnetization case, Ref_s1, in Fig. 7. Here, also, the light curve is shown for a period that is equivalent to ≈20 yr in physical units and the flux values are provided in the units of ergs s⁻¹ cm⁻² Hz⁻¹. In this case, the growth of the instability is not sufficient to show significant structural formation. The shaded regions and the vertical lines indicate different activity states shown by the system, and the spectral behavior at the corresponding times is given in Sect. 4.2.4. A moderately variable emission is noticed in all energy bands except for the R-band and 10¹⁷ Hz emission. Similarly to the Ref_s10 case, a transient feature is observed in the R-band and at 10¹⁷ Hz, due to the localized synchrotron emission (see the red shaded region). However, the energy of the particles responsible for the transient optical emission is less compared to the energy of the particles responsible for the transient X-ray emission. As a result, the optical light curve decays slowly, whereas the X-ray light curve decays rapidly as a consequence of the faster radiative cooling effect.

Fig. 7. Simulated light curve for the Ref_s1 case. The observing frequencies and other details are same as Fig. 6. The flux values from top to bottom are normalized to their maximum values of 2.58 × 10−18, 7.99 × 10−19, 6.65 × 10−23, 1.9 × 10−25, and 3.99 × 10−32, respectively, and given in the units of ergs s−1 cm−2 Hz−1.

4.2.3. Correlation and variability analysis

For correlation studies between these light curves, we performed discrete correlation function (DCF; Edelson & Krolik 1988) for different activity states. This method measures correlation functions without interpolating the data in the temporal domain. The unbinned DCF function is defined as:

$$\begin{array}{c}\hfill {\mathrm{UDCF}}_{\mathit{ij}}=\frac{(a_i-\overline{a})(b_j-\overline{b})}{\sqrt{(({\sigma }_a^2-e_a^2)-({\sigma }_b^2-e_b^2))}},\end{array}$$(11)

where $\overline{a}$ and $\overline{b}$ are the mean of the two discrete data sets a_i and b_j, respectively; σ_a/b and e_a/b are the standard deviation and error associated with each data set. We note that for simplicity, we have not added the simulated errors to the dataset namely, e_a/b is taken as 0 while estimating the UDCF_ij values. Furthermore, the DCF(τ) can be measured by binning the result in time. Averaging over the M pairs for which τ − (δτ)/2 < δt_ij ≪ τ − (δτ)/2, the DCF is given as:

$$\begin{array}{c}\hfill \mathrm{DCF}(\tau )=\frac{1}{M}{\mathrm{UDCF}}_{\mathit{ij}}.\end{array}$$(12)

Figure 8 shows the DCF plots among radio, optical R-band, X-ray (10¹⁷ Hz), and γ-ray (10²⁰ Hz) for the yellow shaded region of the Ref_s10 case, as shown in Fig. 6. We found that the emissions at radio frequencies (1.4 GHz and 43 GHz) show a strong correlation with zero time lag, where each piece of data corresponds to about four months of binned data. This suggests that the emitting regions for the radio bands are co-spatial within that period (Baliyan 2001). Additionally, during this period, the γ-ray emission is correlated with radio emission with 5 t_sc and 2 t_sc lag for 1.4 GHz and 43 GHz, respectively (left panel of Fig. 8). An earlier study by Max-Moerbeck et al. (2014) also found such lag in radio emission from the γ-ray emission. The DCF plots shown in the middle panel of Fig. 8 indicate a correlation between the X-ray emission with radio band emissions with a lag of 3 t_sc and 2 t_sc for 1.4 GHz and 43 GHz, respectively. This suggests that the emission in the X-ray band comes prior to the emission in the radio bands, namely, the highly energetic shocked electrons emit first in the X-ray followed by the radio emission. Similarly, there is a correlation between X-ray and γ-ray with a positive lag of 1 t_sc. Finally, the DCF plots of R-band with all other frequencies are shown in the right panel of Fig. 8. The plots indicate the correlation of R-band emission with radio, X-ray, and γ-ray emission with different time lags. X-ray and R-band emissions exhibit a 1 t_sc lag, whereas the optical emission is 3 t_sc and 2 t_sc ahead of 1.4 GHz and 43 GHz, respectively. The optical and γ-ray emission are strongly correlated with zero lag. We note that during this period of t/t_sc = 20 to 50, rise in the X-ray band emission is followed by the emission at other energy bands. The lags observed between frequencies suggest that the emitting regions are at least cΔτ apart within the ∼4 month period.

Fig. 8. DCF plots between different energy bands for Ref_s10 cases for a time duration t/tsc = 20 to 50 (yellow shaded region in Fig. 6). The vertical lines indicate the lag time when the DCF of each energy band pair peaks.

Similarly, in Fig. 9, the DCF plots are shown for the Ref_s1 case for the red shaded region of Fig. 7. During the X-ray transient phenomena, a strong correlation with zero lag is seen among the radio bands with the γ-ray emission. This suggests that within this duration of t/t_sc = 33 to 45, the radio and γ-ray emitting regions are co-spatial with each other. Similarly correlated emission between radio and γ-ray was previously observed by Ramakrishnan et al. (2015). During the same period, the γ-ray emission is also strongly correlated with R-band emissions with zero lag. However, the optical emission is 1 t/t_sc ahead of radio emissions due to contribution of both synchrotron and EC. In addition, X-ray emission is leading the radio and γ-ray emissions with a positive lag of 4 t/t_sc and 3 t/t_sc, respectively. Further, X-ray emission leads the R-band emission with a time lag of 3 t/t_sc during the X-ray transient activity state. This suggests that the recently shocked particles create another population of highly energetic electrons that emit in X-rays due to the synchrotron process. Afterwards, due to radiative cooling it loose energy and emits in low energy bands. Then, the radiated-out lower energetic population gives the emission in the γ-ray band due to IC scattering.

Fig. 9. DCF plots between different energy bands for Ref_s1 cases for a time duration t/tsc = 33 to 45 (red shaded region in Fig. 7). Vertical lines indicate the lag time when the DCF of each energy band pair peaks.

To quantify the variability in the light curves, we estimated the relative variability amplitude (RVA) or the variability index (Kovalev et al. 2005; Singh et al. 2019), defined as:

$$\begin{array}{c}\hfill \mathrm{RVA}=\frac{F_{\max }-F_{\min }}{F_{\max }+F_{\min }},\end{array}$$(13)

and the uncertainty on RVA is given by:

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{RVA}=\frac{2}{{(F_{\max }+F_{\min })}^2}\sqrt{{(F_{\max }\mathrm{\Delta }{F}_{\min })}^2+{(F_{\min }\mathrm{\Delta }{F}_{\max })}^2},\end{array}$$(14)

where F_max and F_min are the maximum and minimum values of the simulated flux with ΔF_max and ΔF_min uncertainties, respectively. For the purpose of statistical analysis, we have also simulated the error bars for each simulated flux value. It is important to note that these are not the result of simulations or numerical uncertainties. Refer to Acharya et al. (2021) for a further explanation of how these random errors were generated. The estimated RVA value for Ref_s10 and Ref_s1 cases for a time duration of t/t_sc = 20 to 80 at different energy bands is given below in Table 3.

In Ref_s10 case, at radio bands (1.4 GHz and 43 GHz), the variability increases with an increase in the observing frequency. Less variability at lower frequencies can be attributed to the continuous emission coming from the whole plasma column since the emission gets diminished at higher frequencies due to radiative cooling. At higher energy bands, typically in optical and X-ray, localized high synchrotron emission is coming from the shocked and kinked region, resulting in a very high value of RVA. In γ-ray, a gradually decaying nature is observed in the light curve, resulting in a high RVA value. In the lower magnetized case (Ref_s1), at the lower frequencies, the variability amplitude increases. However, in optical and X-ray, a transient phenomenon is observed due to the generation of localized shocks, giving rise to high RVA = 0.99 ± 0.18 and 0.71 ± 0.12 respectively. In the γ-ray band, the variability is similar to the radio bands, owing to a RVA value of 0.51.

Further, from the typical minimum variability timescales (Burbidge et al. 1974), the size of the emitting region in the radio bands is equivalent to the size of the plasma column as the radio emission is coming from throughout the column. However, the minimum variability time-scales obtained for optical and X-ray band is nearly 0.75 t_sc in the Ref_s10 case. This corresponds to the emitting region that covers nearly 30 number of grid cells for both optical and X-ray emission out of (160 × 160 × 240) number of grid cells in the whole computational domain. Similarly, in Ref_s1 case, the size of the emitting region decreases as we progress from 1.4 GHz to 43 GHz. However, the optical and X-ray emissions mostly come from the shocked particles that correspond to the emitting region that encloses approximately 14 and 60 grid cells, respectively. This indicates that such shocked emissions are very localized and mainly come from the kinked or sheared regions.

4.2.4. Spectral energy distribution (SED)

We modeled the multi-wavelength SEDs of the different simulation runs by using the time-dependent fluid and spectral information. Figure 10 provides the time-evolving spectral behavior for the cases Ref_s10 (left) and Ref_s1 (right). The SEDs exhibit two distinct humps at each time, with the low-energy hump caused by the synchrotron and the high-energy hump caused by the EC process.

Fig. 10. Time-evolving broadband spectra for two different cases. Left: simulated SED of Ref_s10 at t/tsc = 26, 33, 39, and 70. Right: simulated SED of Ref_s1 cases at t/tsc = 28, 38, 50, and 70.

For the Ref_s10 case, at the initial time of the instability (i.e., t/t_sc = 26), the synchrotron spectra is flat and the EC spectra shows a gradual decrease in the slope further 10¹⁸ Hz. At this time stamp, an orderly behavior of the magnetic field lines is expected, owing to the slower growth of the perturbation. Therefore, for an observer making a 5° angle with respect to the axis of the column, the toroidal component of the magnetic field would be dominant to provide strong synchrotron emission. In addition, the onset of instability leads to the development of shocks, generating highly energized particles. In addition to the dynamics of the plasma column, such a mechanism can be attributed to the flattening of SED spectra, which results in a broadening of the spread of the individual humps up until higher energies are achieved. At t/t_sc = 33, the shocked particles have cooled down, resulting in marginally lesser emissions at higher frequencies. However, at t/t_sc = 39, few particles have experienced shocks and the light curve shows a comparatively low rise in amplitude. As a result, the spectra becomes flat again at this time stamp. With the evolution of the instability, the magnetic energy dissipates and the total synchrotron flux values at higher frequencies decrease. As a result, the spectra become steep due to the loss of energy at the higher frequencies. At all time stamps, the high energy component of the SED behaves in a similar fashion to the low-energy component, as it has the same electron population responsible for both.

In the right panel of Fig. 10, the multi-wavelength SED for Ref_s1 case is shown for its different activity states, as indicated in Fig. 7. At t/t_sc = 28, a moderately active behavior is observed among the radio bands correlated with γ-ray emission, followed by a transient characteristic feature in the X-ray band at t/t_sc = 38 (see Fig. 7). At this time stamp, the optical emission also shows a rise in the light curve. Such transient features in optical and X-ray bands are an indication of the generation of localized shocks. As the particles pass through shocks, the non-thermal electrons get accelerated to higher energies, depending on the strength of the shocks, and they generate another sub-population of electrons with high electron’s energy. Consequently, the particle spectra and the SED get flattened. Further, due to the high cooling rate, such shocked highly energetic particles cool faster, causing the sharp cut-off in the SED at higher energies. As in the high-magnetization case, at all time stamps, the spectral behavior of both the low- and high-energy components of the SED is similar.

4.2.5. Parametrical comparison of SEDs and Compton dominance

In this work, we have shown a parametrical comparison of SEDs and discussed the effects of different parameters such as σ₀, T and κ on the broad-band spectra. For this purpose, along with σ₀ = 10 and 1, we chose T = 5000 K and 2000 K, κ = 10⁻² and 10⁻³. The detail of each runs is given in Table 2.

In this plasma column setup, the case with a high magnetization value (σ₀ = 10) is more prone to kink instability, whereas in the lower magnetization (σ₀ = 1) case, a throttled growth of kink instability is observed due to the presence of weak magnetic field strength. This results in a comparatively small amplitude of the synchrotron peak. Furthermore, an inadequate dissipation of magnetic energy and a lower extent of loss due to synchrotron prompts the orderly nature of the ensuing magnetic field lines. As a result, a steady emission is observed at the lower radio frequencies up to 10¹³ Hz, with a broad photon spectra, followed by a decrease. Such comparisonal behaviour of SED for σ = 10 and 1 is shown as black (circle marker) and pink (star marker) lines, respectively, in Fig. 11. The Compton dominance (CD) parameter would be one of the intriguing tools that may govern the magnetization of the system. Typically, it is defined as the ratio of the peak of the Compton to the peak of the synchrotron luminosities (Finke 2013; Nalewajko & Gupta 2017). In our work, we calculated the CD in two ways. Firstly, we define CD^peak as the ratio of the Compton peak (high energy hump) to the synchrotron peak (low energy hump) flux densities, noting that the ratio would stay unchanged if we calculated it from the Compton peak to the synchrotron peak luminosities. Furthermore, we defined CD^power as the ratio of EC power to the synchrotron power derived by integrating the flux values across the frequency range. This mainly gives an estimation of how much more EC emission there is than the synchrotron emission. The CD^peak and CD^power values for all the simulation runs are provided in Table 4. For the σ₀ = 10 case, the CD^peak is found to be nearly 1, whereas for the σ₀ = 1 scenario, the CD^peak is estimated to be nearly 10. Such dependencies displayed by Compton dominance on the magnetization parameter have previously been studied by Janiak et al. (2015) and our results are in agreement with their outcomes. Additionally, the CD^power is higher for the lower magnetization case since the EC component of the spectra becomes broader. With change in the temperature of the photon field or the κ values, the energy loss rate of the emitting particles changes. Such effects of photon field temperature and κ values can be seen from the cooling timescale (see Fig. A.1) and the time-evolving particle spectra (refer Fig. A.2). As we decrease the temperature of the photon field, the loss due to EC also decreases. As a result, it demonstrates a substantial effect on the low-energy hump as well as the high-energy hump of the SED. Due to less EC loss, a larger number of highly energetic particles will be present to emit synchrotron emission up to 10¹⁵ Hz. With the considered fluid and spectral information, due to the effects of different radiative losses, the synchrotron emission at higher frequencies get enhanced compared to the higher photon field temperature case. Consequently, the EC emission at the higher end of the spectrum also increases. This results in broadening of the individual component of the spectra. Consequently, even though the CD^peak is found to be ≈0.4, the value of CD^power suggests that the EC power is nearly 2.5 times the synchrotron power. The SED for this low temperature case (Ref_s10_A) is shown as purple dotted (square marker) lines. The blue (plus marker) line represents the SED for the case with a decreased value of κ (Ref_s10_B case) namely, if the distance between the location of the target seed photons and the emitting region increases. Decreasing the κ value affects the loss rate of the emitting particles and, hence, it alters the characteristics of both humps. As in the lower temperature scenario, the synchrotron and EC emission are greater at higher energies, as compared to the high κ values, resulting in a broad photon spectra. The CD^peak for the Ref_s10_B case is estimated to be 0.08 and the CD^power value is about 1.2.

Fig. 11. Parametrical comparison of simulated SEDs for an observer making a 5° angle with respect to the axis of the column estimated at t/tsc = 70.

In summary, it appears that the values obtained from CD^power and CD^peak display a similar trend. However, the estimation of CD^power does not only consider the peak flux density value, but also the shape of the individual spectra that is driven by various radiative cooling processes. As a result, in all cases, the values of CD^power are estimated to be greater than 1, implying greater EC power than synchrotron power.

5. Discussions

The primary focus of this work is to investigate the impact of instability driven shocks on the long-term nature of parsec scale jets in the relativistic regime and their consequences on the multi-wavelength spectra. Previously, Acharya et al. (2021) performed simulations of a helical jet model in the context of long-term variability and studied the dependence of the viewing angle on the observed emission using static particle spectra. In this work, the results pertaining to emission are obtained using time-dependent particle spectra, where the spectral information along with the fluid information evolves with time. It is also important to mention that at parsec-scale lengths, for highly magnetized and relativistic jets, IC scattering due to external photon fields (EC) is evidently more responsible for the observed emission, as compared to IC-CMB. On this basis, we have implemented the EC process in the hybrid framework of the PLUTO code in this work. In addition, our emission modeling approach resembles the multi-zone modeling method, where each grid cell acts as a single emitting region, with fluid spectral information evolving over time.

5.1. Incorporation of EC process

Demonstrating the effect of the EC on the emission signatures of relativistic jets at parsec scale lengths is one of the main objectives of this paper. We note that at this length scale, the effect of IC-CMB emission is significantly lower compared to other high energy emission mechanisms such as external Compton, synchrotron self-Compton, and so on. The parameterizations and approximations used in this work allow for a detailed description of the EC (IC scattering of mono-directional photons) with a ∼1% precision. The executed formalism to obtain emissivities is valid in both the Thomson and Klein-Nishina limits. Our approach includes a Planckian distribution of the target photon field as suggested by Khangulyan et al. (2014) in contrast to δ-function and/or step function approximations for the seed photon fields (Dermer et al. 1992; Dermer & Schlickeiser 1993; Petruk 2009; Finke 2016). In both the Thomson and Klein-Nishina regimes, in the presence of both isotropic and anisotropic radiation fields, the EC mechanism has been rigorously modeled by many authors and those studies include numerically intense calculation (Dermer et al. 1992; Boettcher et al. 1997; Hutter & Spanier 2011; Hunger & Reimer 2016). The analytical calculations performed in Zdziarski & Pjanka (2013) provide a high accuracy level of 0.3% for estimating the interaction rate; however, the usage of special functions in the algorithm inhibits it in its effectiveness and practical usage. First, we characterized the energy losses due to several radiative processes through cooling time-scales (see Fig. A.1). The highly energized electrons lose energy faster and are mainly responsible for the synchrotron emission; whereas, for EC, in the Thomson limit (≈γΘ ≪ 1), electrons lose energy like synchrotron process and in the Klein-Nishina limit (≈γΘ ≫ 1), the loss mechanism follows a logarithmic profile. Furthermore, we have performed 2D toy model simulations that are focused on validating our numerical algorithm and understanding the effects of various parameters on the particle spectra (see Fig. A.2) and multi-band emissions (see Fig. 4).

5.2. Implications of shocks on the emission signatures

To study the multi-wavelength nature of relativistic jets, we simulated a cylindrical plasma column with two different magnetization values and perturbed it with a kink mode instability. The higher magnetization value is more susceptible to the kink mode instability, while in the case of lower magnetization, the throttled growth of kink instability with a mixing of KHI is observed due to the presence of shear and the trans-Alfvénic nature of the flow (Baty & Keppens 2002; Acharya et al. 2021; Acharya & Vaidya 2022). With the evolution of instability, shock formation generates another sub-population of highly energized electrons. In particular, the shocks are formed at the kinked region in the magnetically dominated scenario, and in the lower σ case, shear-driven shocks are produced. Such highly energized electrons are responsible for the observed abrupt rise seen in the multi-wavelength light curves for both cases. A similar remarkable rise in the X-ray light curve has previously been seen in FSRQ PKS 1222+216 (Chatterjee et al. 2021). Very short lived X-ray flares have also been observed in PKS 2005−489 (Zhu et al. 2018), where the variability time-scale is found to be < 30 s. Furthermore, Markowitz et al. (2022) also found modest X-ray flares in Mkn 421 that last for less than a day. The kinematic analysis of high-resolution VLBI images of AGN jets reveals the presence of several blobs that are typically interpreted as recollimation shocks. A combination of multi-frequency observations with 43 GHz VLBI observations suggests that the interaction between travelling shock waves and recollimation shocks could also trigger high energy γ-ray flares (Agudo et al. 2012; Schinzel et al. 2012). Numerical simulations of a conical jet have shown that the interaction between recollimation shock and a traveling shock may produce flaring events in both the single-dish light curves and in the VLBI observations (Fromm et al. 2016).

– One of the important implications of shocks in our study is the appearance of X-ray orphan flares. In both high (Ref_s10) and low (Ref_s1) magnetized cases, a fugacious activity is observed in the X-ray band due to the generation of a sub-population of highly energetic shocked particles. Similar isolated X-ray flares have also been seen before by Abdo et al. (2010) for 3C 273 source due to the presence of localized shocks. During such a transient phase, the X-ray emission comes from very localized regions that lead the emissions from other energy bands. The possible explanation for such a lagged and correlated emission is that the highly energetic shocked particles lose energy due to the radiative cooling mechanisms and emit in the lower energy bands in the presence of tangled magnetic fields via synchrotron process. These lower energetic particles get up-scattered and emit γ-rays via the EC process. In discussing the multi-wavelength correlation study for a large sample of blazars, Liodakis et al. (2018) also offered an analogous explanation. Furthermore, multi-wavelength correlation studies have been carried out by several authors to understand various emission mechanisms and the geometry of emitting locations. For example, Chatterjee et al. (2012), Hovatta et al. (2014), Fuhrmann et al. (2014), Max-Moerbeck et al. (2014), Cohen et al. (2014), Ramakrishnan et al. (2015), and many others have conducted primarily radio, optical, and γ-ray correlation studies with the goal of constraining and locating the size of γ-ray emitting zones. Their finding supports our correlation results in the context of lag and non-lag multi-frequency emission along a timescale of lagged emission, obtained from our simulated light curve, lies between the observed range from several studies.

– Time-evolution modeling of the SED provides a better approach to studying different acceleration mechanisms and constraining the physics of the jet more consistently. In our study, we found that during the transient phase observed in the flux variation, the shocked particles emit in the high frequencies and the broadband spectra get flattened. The flattening of the spectral slope in SED is due to the re-acceleration of electrons, produced as a consequence of localized shocks. Such spectral signatures have also been observed before and studied by Micono et al. (1999), Borse et al. (2021). In addition, such spectral broadening sometimes results in a shift in the peak of the synchrotron component. More recently, a multi-wavelength study of BL Lacertae shows that during the high flux state in X-rays, synchrotron emission is the dominant mechanism for X-ray production when considering that the emitting region is very close to the black hole. Following the X-ray flaring activity, the X-ray component of the SED lies in the higher energy range and occurs due to the EC mechanism (Agarwal et al., in prep.). Similar results were also found before for BL Lacertae by Böttcher et al. (2003).

Particle acceleration due to magnetic reconnection is more effective in studying fast variability scenarios in the case of blazar jets (Narayan & Piran 2012; Giannios 2013; Shukla & Mannheim 2020). However, the current work does not specifically focus on the fast variability nature (Aharonian et al. 2007; Albert et al. 2007; Ghisellini & Tavecchio 2008; Giannios et al. 2009; Barkov et al. 2012; Pryal et al. 2015). Typically, the origin of intra-day and intra-night variable signatures in γ-ray are attributed to magnetic reconnection. Previous studies have demonstrated the role of kink in generating current sheets that could initiate fast magnetic reconnection, providing an efficient way of particle acceleration in AGN and gamma-ray-burst jets (Singh & Mizuno 2016; Bodo et al. 2021; Kadowaki 2021; Medina-Torrejón et al. 2021). However, it is also important to highlight the effects of shocks generated at the shear interface due to the onset of instability on the emission signatures.

5.3. Note on Compton dominance

Magnetic fields threading the central rotating objects are predominantly responsible for the formation of relativistic jets. These jets are converted from initially Poynting flux-dominated flows to matter-dominated flows during their propagation in space (Sikora et al. 2005). Also, MHD instabilities may play a promising role in this process. In this work, we have focused on two different scenarios: one with a dominant kink mode instability and the other one with mixing of both the kink mode and shear-driven KH instability. Both instances differ from each other in terms of the magnetization values. The jet magnetization parameter is not only important for understanding the dynamical evolution of jets but also crucial in governing the dominant particle acceleration mechanism; the Compton dominance parameter can be a measure of the magnetization of the jet (Janiak et al. 2015). Sikora et al. (2005) studied the dependence of jet magnetization on the Compton dominance as a function of the geometry of the emitting region and location of the seed photon source. They found that the Compton dominance is higher for a lower magnetization value, irrespective of the geometry and distance of the photon source from the blazar-emitting region. Similar results on Compton dominance and jet magnetization were also suggested by Nalewajko & Gupta (2017). In our work, from the presented broad-band spectral analysis, we find a Compton dominance, obtained from the peak values of flux densities (CD^peak), of approximately 10 for the σ = 1 case, whereas for the σ = 10 case, it is found to be nearly 1. Such a value of Compton dominance is also observed in FSRQs in contrast to BL Lacs due to the presence of an external photon field. In addition to different magnetization values, we performed simulations with different photon field temperatures and κ values. From our analysis, we found that with the decrease in the temperature of the photon field, the CD^peak decreases and the spread of the spectrum for both humps increases as a result of less EC loss. A similar broadening of individual components of the SED is also observed for all the comparative cases mentioned in Table 4. As a result, the Compton dominance obtained from the power of each components (CD^power) suggests that there is a higher EC power at work compared to the synchrotron power. The CD^peak and CD^power values in the Ref_s10_A case is estimated to be on the order of ≈0.4 and 2.5, respectively. Also, we noticed that the CD^peak and CD^power values decreases with a decrease in the κ value, and it is found to be nearly 0.08 and 1.2, respectively, for the Ref_s10_B case. The values obtained from CD^power and CD^peak seem to follow a similar pattern.

In summary, the synthetic SEDs generated show a presence of mild Compton dominance (CD^peak), particularly for runs with a lower magnetization. The small value can primarily be attributed to the fact that the EC emitting region is irradiated by the photon source that essentially lies beneath (mono-directional photon field). However, if we take into account an isotropic external photon distribution, we may expect to obtain higher CD^peak due to the presence more external photons. Furthermore, isotropic electron distribution is a zeroth-order assumption. According to Kelner et al. (2014), adopting an anisotropic distribution of electrons may have an impact on the measured emission, particularly on the peak of the SED. Implementing such a distribution, however, would be a complex process that is beyond the scope of the current work.

6. Summary

MHD instabilities are one of the vital processes responsible for the observed emission features of blazar jets. The growth of these instabilities lead to the generation of shocks, resulting in the acceleration of particles up to very high energies, which is considered to be one of the prime reasons for the high energy flares seen in these jets. In this work, we performed a parametric study characterizing the EC process and implemented it in the hybrid framework of the PLUTO code. Furthermore, we studied the significance of different radiation mechanisms by performing a numerical simulation of a 2D relativistic slab jet. Lastly, we looked at the multi-wavelength emission of a 3D plasma column in differently magnetized environments where MHD instabilities were present. The key results of this work are summarized as follows.

• We observed a sudden rise in X-ray and optical flux due to the acceleration of emitting electrons near the vicinity of freshly formed shocks generated due to jet instabilities.

• The impact of such localized shocks is manifested in the SEDs, as we find evidence of flattening during the active growth phase of jet instability.

• Subsequent to the passing of shock, the particles cool towards lower energy and their evolution is dominated by losses due to the EC process for the chosen set of parameters. This transformation of dominance in radiating processes gives rise to a change in slope for the SEDs between the optical and X-ray bands.

• The effect of radiative cooling due to the interplay of synchrotron and EC processes is clearly revealed by the shifts in the discrete correlation function. In particular, during the shocked phase, X-ray emission leads all other wave-bands. During the later evolutionary stages, a correlation is observed between the radio and γ-ray emissions owing to the up-scattering of external photons via low energy radio electrons giving rise to γ-ray emission through the EC process.

• The impact of including the EC process is also evident from the observation of the Compton dominance (CD^peak) of ∼10 for the low-magnetization case. Such values for the Compton dominance are typically seen in low-synchrotron peaked (LSP) blazars.

In summary, the present work highlights a unique multi-zone and time-dependent framework for blazar that can qualitatively provide insights on the role played by the interplay of instability-driven shocks in the presence of two radiative loss mechanisms. Furthermore, these jets also show a high degree of linear polarization. Several recent studies have shown that the observed optical polarization variability is correlated with high-energy flaring activity. In a subsequent work, we plan to study polarization signatures of the relativistic jets by performing 3D simulations that are prone to both axisymmetric and non-axisymmetric instabilities.

Acknowledgments

The authors are grateful to Dmitry Khangulyan for the detailed discussion and insightful comments. S.A. is supported by the DST INSPIRE Fellowship and would like to acknowledge the support for her Ph.D. B.V. would like to acknowledge the support from the Max Planck Partner Group Award. I.K.D. thanks Department of Science and Technology’s Swarnajayanti Fellowship (DST/SJF/PSA-03/2016-17) at IISc Bangalore for financial support. The authors are also grateful to the reviewer for his thorough reading and insightful comments, which helped to improve the manuscript. All the computations presented in this work are carried out using the facilities provided at IIT Indore and the Max Planck Institute for Astronomy clusters: VERA, which is a part of the Max Planck Computing and Data Facility (MPCDF).

References

Appendix A: Validation of implementation of EC

In the following, we have discussed the cooling effect of a single emitting particle due to the radiative processes such as synchrotron and external Compton. Furthermore, in Section A.2, particle evolution of a single macro-particle is discussed in the context of the relativistic slab jet problem.

A.1. Calculation of radiative cooling time

Ignoring the loss due to adiabatic expansion, we estimated the cooling time for synchrotron and EC processes for several values of the magnetic field strength, photon field temperature, and κ by separating each of the terms in Equation 1 and treating them individually. It should be noted that while estimating the cooling times for various radiation mechanisms, the loss due to adiabatic expansion was purposely turned off. However, while performing the above-mentioned simulations, we took into consideration the loss owing to adiabatic expansion as well as the loss due to synchrotron, IC-CMB, and external Compton forces. The top panel of Figure A.1 represents the cooling time scales for synchrotron processes only, whereas the bottom two panels represent the cooling time scales for solely EC mechanisms. These time-scales are estimated for an initial power law distribution of particles with energy cutoffs γ_min = 10² and γ_max = 10⁸ as well as a power law index of p = 6. Here, we introduced a parameter, γ_cool, which is defined as an energy for which the radiative cooling timescale is at minimum. The plots of the synchrotron cooling timescale show that the high-energy particles cool faster compared to low-energy particles. Additionally, a stronger magnetic field strength corresponds to stronger cooling and, therefore, the cooling timescale is much smaller. However, γ_cool, remains the same for all B-field values. As seen from the middle and bottom panels, with an increase in T and κ, the cooling time decreases. The γ_cool shifts to a lower energy range with increasing value of T but, on the other hand, it does not have dependence on κ. Such a parametrical behavior of γ_cool is shown in all the panels of Figure A.1 with black dashed lines. The functional form given in Equation 6 implies that at the lower energy range, the loss rate due to EC is similar to that of synchrotron loss. However, in the higher energy range, it has a logarithmic profile (see Equation 34 of Khangulyan et al. 2014). Hence, unlike in a synchrotron, where the maximum energy of the electron is responsible for most of the observed emission or luminosity, in EC, the emission would be maximum when the energy of an electron is equivalent to the energy of the photon. Therefore, the cooling time scale would be minimum for an electron, where the transition occurs from Thomson to the Klein Nishina limit. In other words, the emitting particle that has the minimum cooling time would have the maximum contribution to the observed emission. Furthermore, adopting a particle distribution other than the power-law distribution may alter the cooling timescale profile. However, the typical behavior of the particle in the cooling process would be the same.

Fig. A.1. Cooling time scale of the emitting electron for the synchrotron process only (top panel) and EC process only (middle and bottom panels). The black dotted lines indicate the energy for which the cooling time is minimum. The different line styles correspond to different parameters as stated in the plot.

A.2. Evolution of particle spectra

To calculate the energy loss rate due to several radiative processes, Vaidya et al. (2018) employed an analytical solution to solve Eq. 1 with c₃ = 0. As already mentioned in Section 2, we have included an additional term due to the EC process and solved the Equation 1 numerically by implementing the fourth order Runge-Kutta method (RK4). For validation purposes, we have shown the particle spectra of a single test macro-particle from case 1 of the 2D slab jet problem (see Section 3.1) at the final time of the simulation in the top panel of Figure A.2. In this panel, the cyan solid line depicts the spectra where Equation 1 is solved analytically and the black dashed line represents the spectra where it is solved using the RK4 method. We note that in both the scenarios (without the inclusion of the loss term due to the EC process), the spectra show a great degree of overlap. This validates the correctness of extending our analytical method to the numerical approach. This extension is imperative for including the term related to EC losses.

Fig. A.2. Time evolution of particle spectra for a single test macro-particle in the context of 2D slab jet problem with and without including EC. Top panel shows the spectra at t/tsc = 200 for case 1 with both analytical and numerical solutions.

The energy-loss rate of an electron with and without the inclusion of EC differs significantly from each other. The effect of such energy loss is distinctly visible in Figure A.2, which exhibits the time evolution of particle spectra for a single macro-particle. In our simulations, a single macro-particle refers to a collection of a number of micro-particles or leptons. The second panel of Figure A.2 shows time-evolving particle spectra for a single test macro-particle incorporating loss due to only the synchrotron and IC-CMB mechanisms. The next three panels at the bottom manifest the spectra of the single macro-particle with the incorporation of energy loss due to the EC process, in addition to synchrotron and IC-CMB for different parameters. To understand the dependence of different parameters responsible for the EC process, a set of values for T (= 1000 K and 2000 K) (Donea & Protheroe 2003; Malmrose et al. 2011; Oyabu et al. 2017) and κ (= 10⁻⁷ and 10⁻⁶) were adopted. A noticeable difference is visible between the particle spectra shown in the second and bottom three panels of Figure A.2. The addition of EC loss in Equation 1 results in the cooling of electron in a process that is faster and stronger, compared to the case that includes loss due to synchrotron and IC-CMB only. In addition, higher temperatures and higher values for κ correspond to stronger loss, therefore the γ_max of the emitting particle becomes very small. A similar effect was also seen while considering the behavior of a collection of macro-particles randomly distributed throughout the jet (see Section 3.2).

All Tables

All Figures

	Fig. 1. Cartoon representation of emitting region and external photon fields with the seed photon’s direction. The size and distance of different components are not shown to scale.
In the text

	Fig. 2. 2D representation of particle distribution at t/tsc = 150 (left) and 200 (right) overlaid with the fluid density (grey color). The colour bars show the fluid density and particle unique identities.
In the text

	Fig. 3. Emissivity slices of slab jet for case-1 and case-2 at an observing frequencies ν = R-band and 1017 Hz at t/tsc = 200. Here, the emissivity values are normalized to their individual maximum values.
In the text

	Fig. 4. Emissivity slices of slab jet for case 2 for different observing frequencies at t/tsc = 200. Here, the emissivity values are normalized to their individual maximum values.
In the text

Fig. 5. X − Z cuts of emissivity slices for Ref_s10 case at observing frequencies ν = 1.4 GHz (a), 43 GHz (b), R-band (c), 1017 Hz (d) and 1020 Hz (e) at time stamps t/tsc = 26 (top) and 70 (bottom). The emissivity values are normalized to their maximum values and given in the units of erg s−1 cm−3 Hz−1 str−1. The normalized values of 3.42 × 10−22, 4.23 × 10−22, 3.94 × 10−25, 8.31 × 10−27, 3.4 × 10−35, 2.83 × 10−22, 1.74 × 10−22, 4.76 × 10−26, 6.71 × 10−29, and 3.4 × 10−38 for panels a, b, c, d, e, f, g, h, i, and j, respectively.

In the text

Fig. 6. Simulated light curve for the Ref_s10 case for an observer making 5° angle with the axis of the column. The observing frequencies are ν = 1.4 GHz, 43 GHz, 1.4 × 1014 Hz (R-band), 1017 Hz, and 1020 Hz, respectively, from top to bottom. The vertical black dotted lines and the shaded regions correspond to different activity state of the system (discussed in Sect. 4.2.2). The flux values from top to bottom are normalized to their maximum values 2.64 × 10−17, 2.52 × 10−17, 2.66 × 10−21, 5.13 × 10−24, and 1.54 × 10−31, respectively, and given in units of ergs s−1 cm−2 Hz−1.

In the text

	Fig. 7. Simulated light curve for the Ref_s1 case. The observing frequencies and other details are same as Fig. 6. The flux values from top to bottom are normalized to their maximum values of 2.58 × 10−18, 7.99 × 10−19, 6.65 × 10−23, 1.9 × 10−25, and 3.99 × 10−32, respectively, and given in the units of ergs s−1 cm−2 Hz−1.
In the text

	Fig. 8. DCF plots between different energy bands for Ref_s10 cases for a time duration t/tsc = 20 to 50 (yellow shaded region in Fig. 6). The vertical lines indicate the lag time when the DCF of each energy band pair peaks.
In the text

	Fig. 9. DCF plots between different energy bands for Ref_s1 cases for a time duration t/tsc = 33 to 45 (red shaded region in Fig. 7). Vertical lines indicate the lag time when the DCF of each energy band pair peaks.
In the text

	Fig. 10. Time-evolving broadband spectra for two different cases. Left: simulated SED of Ref_s10 at t/tsc = 26, 33, 39, and 70. Right: simulated SED of Ref_s1 cases at t/tsc = 28, 38, 50, and 70.
In the text

	Fig. 11. Parametrical comparison of simulated SEDs for an observer making a 5° angle with respect to the axis of the column estimated at t/tsc = 70.
In the text

	Fig. A.1. Cooling time scale of the emitting electron for the synchrotron process only (top panel) and EC process only (middle and bottom panels). The black dotted lines indicate the energy for which the cooling time is minimum. The different line styles correspond to different parameters as stated in the plot.
In the text

	Fig. A.2. Time evolution of particle spectra for a single test macro-particle in the context of 2D slab jet problem with and without including EC. Top panel shows the spectra at t/tsc = 200 for case 1 with both analytical and numerical solutions.
In the text