© ESO, 2018

1. Introduction

Extragalactic relativistic jets – collimated outflows of plasma expelled by supermassive black holes residing in central regions of active galactic nuclei and traveling for hundreds of kpc – are among the most fascinating astrophysical structures. Their phenomenology is regulated by complex time-dependent physical processes involving magnetic fields and plasma. Shocks and/or magnetic reconnection sites are able to accelerate charged particles to ultrarelativistic energies, whose presence is flagged by intense non-thermal emission over the entire electromagnetic spectrum, up to the TeV band (e.g., Romero et al. 2017, for a review). Relativistic jets are best studied in blazars (see e.g., Urry & Padovani 1995). The jet of these sources is closely aligned toward the line of sight and, thanks to this favorable geometry, their non-thermal emission is strongly amplified by the relativistic beaming effects.

The observation of the emission from blazar jets at the highest energies, accessible thanks to space (Fermi) and ground-based (Cherenkov arrays) instruments, is revealing numerous unexpected features. One of the most intriguing and challenging aspects is the ultra-fast variability (doubling time ≲ 1 h, down to a few minutes) detected in several blazars (Aharonian et al. 2007; Albert et al. 2007; Arlen et al. 2013; Aleksić et al. 2011, 2014). These very small timescales are often smaller than the light-crossing time of the horizon of the supermassive black hole (e.g., Vovk & Babić 2015) and their interpretation requires extreme physical conditions (e.g., Begelman et al. 2008). In particular, very rapid variability events detected in the most powerful blazars, the flat spectrum radio quasars (FSRQ), are the most challenging to interpret. FSRQ display the classical broad emission lines of quasars, flagging the existence of a broad line region in which clouds of photoionized gas reprocess part of the intense continuum emitted by the accretion flow. Due to the anticipated absorption of the γ-ray radiation above a few tens of GeV through the interaction with soft (UV) photons (e.g., Liu & Bai 2006), one can put robust lower limits on the distance of the emitting region from the central black hole since, to avoid absorption, the emission is constrained to occur beyond the BLR radius R_BLR, of the order of ~ 0.1 pc. For a conical geometry of the jet (with semiaperture θ_j ≈ 0.1) this lower limit on the distance readily translates into a lower limit for the jet radius and (through light-crossing time argument) to the minimum variability timescale, $\mathrm{\Delta }\mathit{t}\gtrsim {\mathit{R}}_{\mathrm{BLR}}{\mathit{\theta }}_{\mathrm{j}}\mathit{/}\mathit{c}\mathrm{𝒟}\mathrm{=}{\mathrm{10}}^{\mathrm{5}}{\mathrm{𝒟}}_{\mathrm{1}}^{-1}$ s (for which we use a Doppler factor¹), a value clearly not compatible with the most extreme minute-timescale events. One of the first FSRQ events showing rapid variability was PKS 1222+122 which, during a flare in June 2010, was detected by the MAGIC telescopes with a flux increasing with a doubling timescale of about 10 min (Aleksić et al. 2011). Since the discovery of PKS 1222+122 other FSRQ have been observed to vary on short timescales. Recently, 3C 279 was detected during a pointed Fermi/LAT observation to vary with a timescale of about 5 min (Ackermann et al. 2016).

In the framework of the classical shock-in-jet scenario for particle acceleration (e.g., Aller et al. 1985) it is clearly hard to explain very short flares. At first sight, the only possibility to explain the rapid variability in FSRQ seems to require the existence of compact (R ≲ 10¹⁵ cm) emission regions embedded in the jet at distances of the order of ~0.1 pc from the central engine (e.g., Ghisellini & Tavecchio 2008; Tavecchio et al. 2011). Generally, such compact regions have been identified in active substructures of the jet, such as plasmoids resulting from efficient reconnection of the magnetic field (Giannios 2013; Petropoulou et al. 2016) or in turbulent cells embedded in the relativistic flow (Narayan & Piran 2012; Marscher 2014). These scenarios have been widely discussed in the past (see also Aharonian et al. 2017). However, one can adopt a different view and explore the possibility that the emission occurs within the entire (or a large fraction of the) jet. This could be reconciled with the fast variability if the jet were to suffer a strong recollimation interacting with an external medium. Jet recollimation has been widely studied in the past (e.g., Komissarov & Falle 1998; Sokolov et al. 2004; Stawarz et al. 2006; Nalewajko & Sikora 2009). Bromberg & Levinson (2009) extended the calculations including the effect of strong radiative losses of the shocked plasma. The cooling, causing a loss of pressure in the compressed region, allows the jet to be “squeezed” by large factors. Bromberg & Levinson applied their scenario to the case of BL Lac objects and radiogalaxies but – as proposed in Tavecchio et al. (2011) – it can also be applied to FSRQ, for which the cooling of the plasma (dominated by the inverse Compton scattering of the ambient soft photons (e.g., Ghisellini et al. 1998) is expected to be even more severe than for other sources.

In this paper we intend to study the recollimation of a jet subject to important radiative losses extending the previous approximated semi-analytical approach of Bromberg & Levinson (2009), by performing high-resolution axisymmetric numerical simulations using the state-of-the-art numerical code PLUTO (Mignone et al. 2007, 2012). The structure of the paper is as follows. In Sect. 2 we illustrate the model and the setup used, in Sect. 3 we show the results of simulations for stationary states and for the time-dependent evolution of a density perturbation, and in Sect. 4 we conclude.

2. The model

We want to construct a numerical model of the confinement of a relativistic jet by the pressure of the surrounding medium. We assume axisymmetry and make use of the cylindrical coordinates r and z, where z is the coordinate along the jet axis and r is the radial coordinate. We inject a conical relativistic jet, with opening angle θ₀, in an overpressured confining region with a pressure profile ${\mathit{p}}_{\mathrm{ext}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathit{p}}_{\mathit{a}}{\left(\frac{\mathit{z}}{{\mathit{z}}_{\mathrm{0}}}\right)}^{-2}\mathit{,}$(1)and a density profile ${\mathit{\rho }}_{\mathrm{ext}}\mathrm{=}{\mathit{\rho }}_{\mathit{a}}{\left(\frac{\mathit{z}}{{\mathit{z}}_{\mathrm{0}}}\right)}^{-2}\mathit{,}$(2)where z₀ represents the height at which the jet starts to interact with the confining region.

The jet is injected in the confining region (z₀ = 1) with a Lorentz factor γ_j and a radius r_j = θ₀z₀, and is subject to radiative losses by non-thermal processes. In this paper we do not consider the dynamical effects of a magnetic field, therefore the evolution equations determining the jet dynamics are the continuity equation ${\mathit{\partial }}_{\mathit{\nu }}\mathrm{\left(}\mathit{\rho }{{u}}^{\mathit{\nu }}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}$(3)where ρ is the proper density and u^μ is the four-velocity, and the energy momentum conservation ${\mathit{\partial }}_{\mathit{\nu }}{\mathit{T}}_{\mathit{\mu \nu }}^{}\mathrm{=}{\mathit{S}}^{\mathit{\mu }}\mathit{,}$(4)where ${\mathit{T}}_{\mathit{\mu \nu }}^{}$ is the stress-energy tensor of a perfect gas ${\mathit{T}}_{\mathit{\mu \nu }}^{}\mathrm{=}\mathit{w}{\mathit{u}}^{\mathit{\mu }}{\mathit{u}}^{\mathit{\nu }}\mathrm{+}\mathit{p}{\mathit{g}}^{\mathit{\mu \nu }}\mathit{,}$(5)where w and p denote, respectively, the proper enthalpy and pressure, and g^μν is the metric tensor for a flat space. S^μ denotes a source term associated with radiative losses, whose form is specified below. The system of Eqs. (3) and (4) is completed by providing an equation of state relating w, ρ, and p. Following Mignone et al. (2005), we adopt the following prescription: $\mathit{w}\mathrm{=}\frac{\mathrm{5}}{\mathrm{2}}\mathit{p}\mathrm{+}\sqrt{\frac{\mathrm{9}}{\mathrm{4}}{\mathit{p}}^{\mathrm{2}}\mathrm{+}{\mathit{\rho }}^{\mathrm{2}}}\mathit{,}$(6)which closely reproduces the thermodynamics of the Synge gas for a single-species, relativistic perfect fluid with a smooth transition from the adiabatic exponent Γ = 5/3 in the non-relativistic limit to Γ = 4/3 in the ultra relativistic case (here and in the following we always put c = 1).

For radiative losses we assume a very simple form; they are taken to be proportional to the pressure and we additionally assume that it is only the gas above a certain threshold temperature T_c that radiates. The first assumption is widely adopted (e.g., Komissarov & Falle 1998; Bromberg & Levinson 2009). For the case under study here the cooling is dominated by the inverse Compton emission on a fixed target radiation field. The emissivity (the frequency integrated radiated power per unit volume) of the non-thermal particles measured in the plasma comoving frame, ϵ ≡ −S⁰, can be expressed as: $\mathit{ϵ}\mathrm{=}\frac{\mathrm{4}}{\mathrm{3}}{\mathit{\sigma }}_{\mathrm{T}}\mathit{c}\hspace{0.17em}{\mathit{U}}_{\mathrm{ext}}{\mathrm{\int }}_{{\mathit{\gamma }}_{\min }}^{{\mathit{\gamma }}_{\max }}\mathit{n}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}{\mathit{\gamma }}^{\mathrm{2}}\mathrm{d}\mathit{\gamma ,}$(7)where σ_T is the Thomson cross section, U_ext is the energy density of the external radiation field, γ is the electron Lorentz factor and n(γ) is the non-thermal electron energy distribution, for simplicity assumed to be a power law with slope 2, n(γ) = kγ^-2, in the range γ_min<γ<γ_max, with γ_min ≪ γ_max (all quantities are expressed in the flow frame). Performing the integral one obtains ϵ = (4/3)σ_TcU_extkγ_max. The energy density of the non-thermal relativistic electrons is: ${\mathit{U}}_{\mathrm{e}}\mathrm{=}{\mathit{m}}_{\mathit{e}}{\mathit{c}}^{\mathrm{2}}{\mathrm{\int }}_{{\mathit{\gamma }}_{\min }}^{{\mathit{\gamma }}_{\max }}\mathit{n}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}\mathit{\gamma }\mathrm{d}\mathit{\gamma }\mathrm{=}\mathit{k}{\mathit{m}}_{\mathit{e}}{\mathit{c}}^{\mathrm{2}}\mathrm{\Lambda }\mathit{,}$(8)where Λ = ln(γ_max/γ_min). Combining the two equations above, we can write: $\mathit{ϵ}\mathrm{=}\frac{\mathrm{4}{\mathit{\sigma }}_{\mathrm{T}}}{\mathrm{3}{\mathit{m}}_{\mathit{e}}\mathit{c}}\frac{{\mathit{\gamma }}_{\max }}{\mathrm{\Lambda }}\hspace{0.17em}{\mathit{U}}_{\mathrm{ext}}{\mathit{U}}_{\mathrm{e}}\mathrm{=}\frac{{\mathit{p}}_{\mathrm{e}}}{{\mathit{\tau }}_{\mathrm{c}}}\mathit{,}$(9)where in the last step we use p_e = U_e/ 3 (valid for ultrarelativistic electrons) and we have defined an effective cooling time: ${\mathit{\tau }}_{\mathrm{c}}\mathrm{=}\frac{{\mathit{m}}_{\mathit{e}}\mathit{c}\hspace{0.17em}\mathrm{\Lambda }}{\mathrm{4}{\mathit{\sigma }}_{\mathrm{T}}{\mathit{\gamma }}_{\max }{\mathit{U}}_{\mathrm{ext}}}\mathrm{·}$(10)The pressure of the non-thermal particles is supposed to be a fraction ξ_e of the thermal one, that is, p_e = ξ_ep. We note that the relation ϵ ∝ p_e is valid if, as assumed here, the electron losses are determined by the inverse Compton emission on a fixed target radiation field. In case where the losses are dominated by synchrotron or synchrotron-self Compton emission one should also consider the role of the magnetic field. The assumption on the critical temperature is intended to mimic the idea that high temperatures flag the presence of relativistic particles in the flow. A final point concerns the fact that (as also pointed out by Nalewajko & Sikora 2009) in a realistic treatment one should assume that relativistic non-thermal particles are accelerated close to the shock fronts and then advected in the other regions of the jet. Both the electron distribution and ξ_e should thus be treated as functions of the position in the flow. For simplicity this is not considered here. In any case, we expect that the global dynamical effects of cooling, such as those explored here, do not dramatically depend on these details. Similarly to Bromberg & Levinson (2009) the quantities characterizing the non-thermal population are assumed constant in the jet. Implicitly this choice assumes the existence of some mechanism (e.g., turbulence) continuously supplying the energy lost by particles through the emission.

With the assumptions above, and assuming that, in the rest frame of the gas, radiation is isotropic (for our applications we can safely neglect the anisotropy of the external Compton radiation field, Dermer (1995), Ghisellini & Tavecchio (2010)), we can express S^μ as ${\mathit{S}}^{\mathit{\mu }}\mathrm{=}\left(\mathrm{-}\frac{{\mathit{\xi }}_{\mathrm{e}}\mathit{p}}{{\mathit{\tau }}_{\mathrm{c}}}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\right)\mathit{.}$(11)In the laboratory frame we can thus write ${\mathit{S}}^{\mathit{\mu }}\mathrm{=}\mathrm{-}\frac{{\mathit{\xi }}_{\mathrm{e}}\mathit{p}}{{\mathit{\tau }}_{\mathrm{c}}}{\mathit{u}}^{\mathit{\mu }}\mathit{.}$(12)Equations (3) and (4) are solved numerically by using the adaptive mesh refinenment (AMR) version of the PLUTO code (Mignone et al. 2007, 2012), with a second order scheme and HLLC Riemann solver (Mignone & Bodo 2005). We perform two-dimensional axisymmetric simulations using the cylindrical coordinates r and z on a grid that covers the domain 0 <r<r_D, 1 <z<z_D, where r_D = 1.5z₀ and z_D = 9z₀. We make use of six levels of refinement, the base grid is made up of 64 × 384 points and we have an equivalent maximum resolution of 4096 × 24 576 points. The boundary conditions are reflective on the axis r = 0, outflow at the outer boundaries z_D = 9z₀ and r_D = 1.5z₀. At the lower boundary, z₀ = 1, and for r<r_j we have inflow conditions injecting the jet flow, while, for r>r_j, we have reflective conditions.

The jet is injected with a proper density ρ_j = 10^-6ρ_a, where ρ_a is the ambient density at z = z₀, with a Lorentz factor γ_j = 10 and pressure is p_j = 7 × 10^-3p_a. The corresponding energy flux is: ${\mathit{L}}_{\mathit{j}}\mathrm{=}\mathit{\pi }{\mathit{\theta }}_{\mathrm{0}}^{\mathrm{2}}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{w}}_{\mathit{j}}{\mathit{\gamma }}_{\mathit{j}}^{\mathrm{2}}{\mathit{v}}_{\mathit{j}}\mathit{,}$(13)where w_j, according to Eq. (6), is ${\mathit{w}}_{\mathit{j}}\mathrm{=}\frac{\mathrm{5}}{\mathrm{2}}{\mathit{p}}_{\mathit{j}}\mathrm{+}\sqrt{\frac{\mathrm{9}}{\mathrm{4}}{\mathit{p}}_{\mathit{j}}^{\mathrm{2}}\mathrm{+}{\mathit{\rho }}_{\mathit{j}}^{\mathrm{2}}}\mathit{.}$(14)In the simulations we fix the parameters to values representative for the flaring state of PKS 1222+216 (e.g., Tavecchio et al. 2011). For the energy flux we use L_j = 10⁴⁶ erg s^-1. Using γ_max = 10⁵, Λ ≃ 10 (quite insensitive on the exact value of γ_min) and U_ext = 3 × 10^-2 (suitable to reproduce the IR radiation field of the torus (e.g., Ghisellini & Tavecchio 2008) we obtain τ_c ≃ 3 × 10⁴ s. Since we assume that the jet is already beyond the BLR, we fix z₀ to the BLR radius, z₀ = 7 × 10¹⁷ cm. The cooling time is thus of the order of τ_c = 10^-3z₀/c. We investigate the effects of radiative losses by performing a series of simulations with different values of ξ_e in the range 0.01−0.1.

3. Results

3.1. The equilibrium configuration

The characteristics of recollimation shocks resulting from the interaction of an under-pressured relativistic jet with external confining material have been discussed in the past; in particular by Komissarov & Falle (1998), Bromberg & Levinson (2007, 2009), Nalewajko & Sikora (2009). The gross structure consists of a contact discontinuity separating the shocked jet layer and the ambient medium and a recollimation shock. Komissarov & Falle (1997) derived, under suitable approximations, simple analytical formulae for the geometrical properties of the shock. Assuming a power law profile for the pressure of the external gas, p_ext(z) = az^{− η}, the profile r(z) of the reconfinement shock follows the differential equation: $\frac{\mathrm{d}\mathit{r}}{\mathrm{d}\mathit{z}}\mathrm{=}\frac{\mathit{r}}{\mathit{z}}\mathrm{-}\mathit{A}{\mathit{z}}^{\mathit{\delta }}\mathit{,}$(15)where δ = 1−η/ 2 and $\mathit{A}\mathrm{=}\mathrm{(}\mathit{\pi a}{\mathit{\theta }}_{\mathrm{0}}^{\mathrm{2}}\mathit{c}\mathit{/}\mathit{\mu }{\mathit{\beta }}_{\mathit{j}}{\mathit{L}}_{\mathit{j}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$, where μ = 0.7. The solution of this equation, with the condition r = r₀ at z = z₀ is: $\mathit{r}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathit{z}{\mathit{\theta }}_{\mathit{j}}\mathrm{-}\frac{\mathit{A}}{\mathit{\delta }}\mathit{z}{\mathit{\theta }}_{\mathit{j}}\mathrm{[}{\mathit{z}}^{\mathit{\delta }}\mathrm{-}{{\mathit{z}}_{\mathrm{0}}^{\mathit{\delta }}}^{\mathrm{]}}\mathit{,}$(16)for η ≠ 2, and r(z) = zAln(z_c/z), for η = 2 (i.e., δ = 0) where z_c = z₀exp(r₀/z₀A). For the contact discontinuity surface a similar analytical expression can be derived (Bromberg & Levinson 2007), r_c(z) = r₀(z/z₀)^{η/ 4}.

These analytical solutions do not include the possible cooling of the post-shock plasma. Bromberg & Levinson (2007) and Bromberg & Levinson (2009) constructed a class of semianalytical models for the confinment of the jet by the external pressure including also the effect of the pressure loss due to radiation losses. As above, these solutions are subject to some approximations and limitations: in the shocked layer the flow parameters are assumed to depend only on the coordinate along the jet axis z (gradients along r have been included by Nalewajko & Sikora 2009) and, more critically, they do not get the structure of the solutions for values of z>z^∗, where z^∗ is the point where the recollimation shock reaches the jet axis.

Here we use a different approach, by which we overcome the two limitations mentioned above, solving numerically the time dependent equations and following the jet evolution until it reaches a steady configuration. We start the simulations at z = z₀ with a conical jet of opening angle θ₀ established over the entire computational domain. As the simulation evolves in time, the jet is compressed by the higher pressure of the ambient medium, the recollimation shock is formed and the system evolves towards its steady structure, that is reached at about t = 300z₀/c. The main features of the steady solutions mentioned above can be observed in Fig. 1, where we show the Lorentz factor distribution for the case with the smaller value of ξ_e = 10^-2. In the figure we see the region of unshocked jet material separated by the recollimation shock from the shocked jet layer. A contact discontinuity separates this shocked material from the external medium, which appears in black. The conical recollimation shock reaches the jet axis at z = z^∗ ~ 3z₀, where it is reflected and gives rise to a strong jet deceleration. From the detailed view of the pressure distribution, shown in the left panel of Fig. 2, we can observe the formation of a complex shock structure, which leads to a pressure maximum for 3.6z₀<z< 3.8z₀. Having linked the emissivity of non-thermal particles to the value of the pressure (Eq. (11)), the region of maximum pressure corresponds also to a maximum of radiative losses. Therefore, in the region around the axis, the jet, having lost a fraction of its energy flux by radiation, will remain at a Lorentz factor quite a lot lower than its injection value (Fig. 1). This slower region on the axis is surrounded by a faster region, formed by material that has suffered less radiative losses.

Considering a jet propagating towards the observer with viewing angle θ_v = 1 /γ_j, the Doppler factor will be and, with the assumptions described above, the apparent (Doppler boosted) emissivity of each fluid element as measured by the observer is given by: ${\mathit{ϵ}}_{\mathrm{obs}}\mathrm{=}\frac{{\mathit{p}}_{\mathrm{e}}}{{\mathit{\tau }}_{\mathrm{c}}}\hspace{0.17em}{\mathrm{𝒟}}^{\mathrm{3}}\mathit{.}$(17)The location of the emissivity maximum depends both on the distribution of pressure and on the distribution of Lorentz factor within the flow. Since the Lorentz factor is lower on the axis and increases away from it, we expect to find the maximum of ϵ_obs away from the axis. In fact from the right panel of Fig. 2, where we display the distribution of the observed emissivity, we can see that most of the observed luminosity originates from an elongated region at some distance from the axis, where the product of the pressure and the Lorentz factor is maximized. We remark that, although the thickness of this region is very small (< 0.01 z₀), its length is of the order of ≈ 0.2 z₀.

Fig. 1 Distribution of the Lorentz factor for the case with ξe = 0.01, showing the main features of the steady solution. The red rectangle individuates the region shown in Fig. 2. We notice that the radial scale is strongly stretched.

Fig. 2 Left panel: pressure distribution in the region marked by the red rectangle in Fig. 1 for the case with ξe = 0.01. The arrows indicate the main features already displayed in Fig. 1. Right panel: distribution of the observed emissivity (see Eq. (17)

Fig. 3 Distribution of the Lorentz factors for the three cases with different values of ξe, that is ξe = 0.01 (left panel), ξe = 0.03 (middle panel), and ξe = 0.05 (right panel).

The effects of increasing the radiative losses are illustrated by Fig. 3, where we show the distributions of the Lorentz factor for three cases with different values of ξ_e. From the Figure we can observe a progressive strong reduction in transverse size of the region of shocked jet material between the recollimation shock and the contact discontinuity, for increasing ξ_e. The “squeezing” of the jet also becomes more effective and the position of the minimum jet transverse radius comes closer to the point where the recollimation shock reaches the axis. The Lorentz factor on the axis is also more strongly decreased because of the increased energy losses. We note also that after a moderate re-expansion, following the development of the reflection shock, the jet is recollimated again.

Figure 4 shows a zoomed view of the distributions of pressure and observed emissivity for the case with ξ_e = 0.05. The pressure distribution shows an overall shape similar to the previous case, shown in Fig. 2, however it is much reduced in size and the maximum value is about one order of magnitude larger. Consequently, the emission also has a similar shape, with much reduced size and higher peak value.

In order to make this result more quantitative, in Fig. 5 we plot the radial width w of the emission region at the point where the value is half the maximum as a function of ξ_e. We observe a progressive decrease of the width for increasing values of the energy transferred to the non-thermal component. For the largest ξ_e = 0.05 the width is of the order of 10^-3z₀, slightly smaller than the estimate of Bromberg & Levinson (2009). Therefore, even in non-extreme conditions, the focusing due to the recollimation is so effective that the brightest region of the jet is as small as ≈ 7 × 10¹⁴ cm, not far from the size estimated for the rapid variability observed in 1222+216 (Aleksić et al. 2011; Tavecchio et al. 2011).

Fig. 4 Left panel: pressure distribution in the region marked by the red rectangle in Fig. 1 for the case with ξe = 0.05. Right panel: distribution of the emissivity in the observer frame (see Eq. (17).

For illustration, we also report in Fig. 6 the profiles (integrated in planes normal to the jet axis) along the jet axis of the observed emissivity (left panel) $\mathit{L}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{\pi }{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{r}}_{\mathit{D}}}{\mathit{ϵ}}_{\mathrm{obs}}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathit{r}\mathrm{d}\mathit{r,}$(18)where ϵ_obs is defined in Eq. (17) and of the jet energy flux (right panel) $\mathit{P}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{\pi }{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{r}}_{\mathit{D}}}\mathit{w}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}{\mathit{\gamma }}^{\mathrm{2}}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}{\mathit{v}}_{\mathit{z}}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathit{r}\mathrm{d}\mathit{r,}$(19)for the three values of ξ_e. Both quantities are normalized to the jet power. The overall shape of the integrated emissivity L(z) (left panel) is the same for all values of the losses, displaying a broad maximum around z₀ ~ 1.5 in correspondence with the onset of the recollimation shock and a second, narrower maximum marking the region of minimum transverse size. For the lowest value of ξ_e = 10^-2 (black line) this second peak is relatively broad and not very high. For larger values of ξ_e (red and green lines) the overall emissivity increases and the second maximum becomes quite narrow, flagging the compactness of the emission region. We note that, in all cases, the area under the two peaks is of the same order, indicating that the total received luminosity comprises similar contribution from the large (hence slowly varying) recollimation shock and the compact (rapidly variable) region. The profiles of the energy flux (right panel) show the increasing importance of energy losses for increasing values of ξ_e. After a smooth decrease, marking the progressive development of the recollimation shock, the curves (especially the two with ξ_e = 0.03 and 0.05) display a sudden jump, corresponding to the increased losses close to the region of smaller radius. After the expansion produced by the reflected shock, the energy flux remains approximately constant. The fraction of the jet energy flux lost through the emitted radiation goes from 15% for ξ_e = 0.01 to a fraction of about 30% for ξ_e = 0.05. These values are in agreement with the results of BL09.

Fig. 5 Transverse size of the emission region as a function of ξe.

Fig. 6 Left panel: observed emissivity integrated over planes normal to the z axis as a function of the z coordinate. Right panel: energy flux integrated over planes normal to the z axis as a function of the z coordinate.

It is also interesting to consider the efficiency ϵ_diss of the dissipation of the kinetic flux of the jet that we define, similarly to Nalewajko & Sikora (2009), as L_kin = ^∫ρc²v_zγ_j(γ_j−1), where the integral is performed on planes normal to the z direction. Comparing the values of L_kin at injection and after the recollimation region, we derive ϵ_diss = 18% for ξ_e = 0.01 and ϵ_diss = 35% for ξ_e = 0.05. These values are larger than the average shock dissipation efficiency derived by Nalewajko & Sikora (2009), that found values lower than 10% for γ_jθ_j = 1. The reason for this higher dissipation efficiency should be looked for in the important radiation losses of the system, not included in the Nalewajko & Sikora treatment. In fact, radiation losses lead to a lower value of the gas pressure downstream from the recollimation shock. In order to maintain pressure balance with the external material, the system has thus to increase the fraction of jet luminosity dissipated into pressure. This is indeed naturally achieved through a narrower profile of the recollimation shock which increases the dissipation because of the larger component of the velocity normal to the shock.

3.2. Evolution of perturbations and variability

In the previous section we studied in detail the stationary equilibrium structure of the recollimation shock. For astrophysical applications, it is important to have a glimpse on the temporal behavior of the structure and its properties when the jet is subject to perturbations.

Fig. 7 Pressure distribution (logarithmic scale) for the simulation including a propagating density perturbation. The three panels show the pressure in the stationary (equilibrium) state, after a time t = 0.9z0/c and, on the left, at t = 5z0/c. The figure clearly shows the high-pressure wake induced by the perturbation in the jet. In the last panel, the inward shift of the recollimation nozzle is clearly visible.

Fig. 8 Profiles (integrated in planes normal to the jet axis) along the jet axis of the observed emissivity calculated for different times during the propagation of the perturbation. In both panels the red lines display the unperturbed profile. The left panel shows the profile a short time after the injection, while in the right panel we report the profiles at later times.

In order to study the variability behavior, we injected a perturbation at the jet inlet and followed its evolution as propagated along the jet. We choose to change the jet density, keeping a constant temperature, with a perturbation of the form: $\mathit{\rho }\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\left\{\mathrm{1}\mathrm{+}\mathit{ϵ}\hspace{0.17em}\exp \left[\mathrm{-}{\left(\frac{\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}}{\mathit{\tau }}\right)}^{\mathrm{2}}\right]\right\}\mathit{,}$(20)where ϵ is the relative amplitude of the perturbation, τ is the temporal width of the perturbation and t₀ is the time at which the perturbation reaches its maximum value. From now on we take t₀ as the origin of time. For τ we took the value 0.05z₀/c and we explored different values of ϵ. Our results show that only quite large values of ϵ give appreciable effects on the observed emissivity; for example ϵ = 1 gives only a 10% increase in the emissivity, therefore we discuss the results obtained for ϵ = 10. In Fig. 7 we show the logarithmic pressure distributions at three different times. In the left panel we have the equilibrium distribution. The middle panel is for t = 0.9z₀/c and we can see the propagating perturbation that leaves, on its sides, a high pressure wake that compresses and modifies the structure of the recollimation shock. Finally the right panel is for t = 5z₀/c, the perturbation is outside the region displayed in the Figure, however, its wake is still visible and the recollimation shock is still deformed; the point at which the recollimation shock reaches the axis has moved inward.

In Fig. 8 we show the profiles (integrated in planes normal to the jet axis) along the jet axis of the observed emissivity (Eq. (17)) for different times. The left panel shows the first phases of propagation (corresponding to t = 0.21,0.67and1.35z₀/c), with the red curve showing the unperturbed profile for comparison. The peak corresponding to the perturbation decreases its amplitude as it moves outward. This decrease is related to the jet expansion. Because of the higher pressure and the shock deformation in the wake of the perturbation (see Fig. 7), the inner jet maintains an emissivity that is approximately twice the unperturbed value even at larger times. In the right panel, the red curve always shows the unperturbed distribution for comparison, while the blue curve, corresponding to the time at which the perturbation reaches the recollimation region (t = 1.8z₀/c), shows a burst of emissivity, with an observed luminosity increasing by a factor of ≈ 5 with respect to the stationary case. Finally, the green curve in the right panel corresponds to a time when the perturbation has moved outside the displayed region, strongly decreasing its amplitude. In this case the peak in the emissivity at the recollimation point has moved inward, following the shift of the recollimation nozzle visible in the maps in Fig. 7.

These results, although obtained under rather idealized assumptions, outline a coherent scenario, qualitatively compatible with the behavior of powerful blazars during active states. In fact, the light curves of FSRQ during typical flare episodes show a long-term (weeks) increase of the flux (that in our simulations would be provided by the long duration active state of the perturbed inner jet), on top of which short-timescale (days) variations are often detected. These rapid variations could be related to the excitation of the recollimation nozzle by the incoming perturbation. In real astrophysical conditions we could also imagine a series of perturbations with different properties injected into the flow that would account for the erratic behavior of the real light curves.

4. Discussion and conclusions

We have developed a series of axisymmetric simulations of a relativistic jet subject to recollimation caused by an external medium using the state-of-the-art code PLUTO and including the radiative losses mediated by a population of non-thermal electrons in the jet. The set up of the simulations is tailored to model jets associated to FSRQ, in particular the source PKS 1222+216, displaying rapid variability at VHE. Our simulations overcome several of the approximations introduced in the semi-analytical approach of Bromberg & Levinson (2009).

The radiative cooling is then assumed to be dominated by the inverse Compton scattering of external radiation (Sikora et al. 1994; Ghisellini & Tavecchio 2009). The estimated effective timescale for radiative losses, τ_c, is much lower then the dynamical timescale τ_dyn ~ z_o/c, ensuring the effective cooling of the shocked plasma. In order to allow a sizeable energy loss by the jet we have also assumed that the emitting electrons, while cooling, are continuously heated, efficiently converting the energy dissipated by the shock into radiation. We have shown that, even with a relatively minor contribution of the non-thermal particles to the total pressure of the shocked material, p_e/p = ξ_e = 0.05, the fraction of energy flux lost by the flow can be as large as 30%, resulting in a significant deceleration of the flow downstream the recollimation zone.

The significant pressure loss caused by the radiative cooling allows the jet to be focused and extremely “squeezed” by the external pressure. Consistently with the early findings by Bromberg & Levinson (2009), our simulations show that the minimum radius of the recollimation region can be as small as 10^-3z₀, implying a light crossing time comparable to the timescale observed for PKS 1222+216. We have also produced maps for the observed emissivity of the flow, including the effects of the Doppler beaming. These show that the emission is concentrated in a very thin elongated region close to the jet axis. We need to caution that the release of the axisymmetry constraint may lead to instabilities of this narrow region, whose outcome can however be assessed only through three-dimensional simulations.

We have also explored the behavior of the recollimation shock in response to density perturbations injected at the base of the jet. Although highly idealized, these simulations show a rather complex phenomenology, that could be, at least qualitatively, related to the typical flaring pattern of powerful FSRQ.

We would like to mention that our calculations, similarly to those of BL09, should be considered somewhat optimistic. In fact the radiative losses induced in the system are maximized by the specific form of the cooling term adopted in the simulations, simply proportional to the thermal gas pressure. As already noted in Sect. 2, such prescription implicitly assumes that the emitting electrons are continuously heated by some unspecified mechanism, allowing to radiate an amount of energy larger than that put into the non-thermal component just after the shock front. As discussed by BL09, possible heating mechanisms include non-linear oscillations of the contact discontinuity surface or acceleration by shear flows. The investigation of these crucial aspects could be the subject of future specific simulations.

Acknowledgments

This work has been partly funded by a PRIN-INAF 2014 grant, the PRIN INAF SKA-CTA grant and a PRIN MIUR 2015. We acknowledge the CINECA award under the ISCRA initiative, for the availability of high-performance computing resources.

References

All Figures

	Fig. 1 Distribution of the Lorentz factor for the case with ξe = 0.01, showing the main features of the steady solution. The red rectangle individuates the region shown in Fig. 2. We notice that the radial scale is strongly stretched.
In the text

	Fig. 2 Left panel: pressure distribution in the region marked by the red rectangle in Fig. 1 for the case with ξe = 0.01. The arrows indicate the main features already displayed in Fig. 1. Right panel: distribution of the observed emissivity (see Eq. (17)
In the text

	Fig. 3 Distribution of the Lorentz factors for the three cases with different values of ξe, that is ξe = 0.01 (left panel), ξe = 0.03 (middle panel), and ξe = 0.05 (right panel).
In the text

	Fig. 4 Left panel: pressure distribution in the region marked by the red rectangle in Fig. 1 for the case with ξe = 0.05. Right panel: distribution of the emissivity in the observer frame (see Eq. (17).
In the text

	Fig. 5 Transverse size of the emission region as a function of ξe.
In the text

	Fig. 6 Left panel: observed emissivity integrated over planes normal to the z axis as a function of the z coordinate. Right panel: energy flux integrated over planes normal to the z axis as a function of the z coordinate.
In the text

Fig. 7 Pressure distribution (logarithmic scale) for the simulation including a propagating density perturbation. The three panels show the pressure in the stationary (equilibrium) state, after a time t = 0.9z0/c and, on the left, at t = 5z0/c. The figure clearly shows the high-pressure wake induced by the perturbation in the jet. In the last panel, the inward shift of the recollimation nozzle is clearly visible.

In the text

	Fig. 8 Profiles (integrated in planes normal to the jet axis) along the jet axis of the observed emissivity calculated for different times during the propagation of the perturbation. In both panels the red lines display the unperturbed profile. The left panel shows the profile a short time after the injection, while in the right panel we report the profiles at later times.
In the text