© G. Fichet de Clairfontaine et al. 2021

1. Introduction

Emission from relativistic jets in radio-loud active galactic nuclei (AGN) has been detected from the radio band up to the teraelectronvolt range and shows frequent episodes of high flux variability in many sources. The emission is generally ascribed to a population of relativistic particles inside the jet that interact with magnetic and photon fields to produce nonthermal radiation over a large wavelength range. In blazar-type AGN, the emission from the jet plasma is amplified in the observer frame by relativistic beaming effects as the jet axis is aligned with the line of sight, while the observed emission is weaker in radio galaxies, where the jet is seen under a larger angle. The multiwavelength emission of AGN requires a mechanism able to reaccelerate particles as they travel along the jet (Blandford & Koenigl 1979). The processes that are most frequently proposed are acceleration from internal shocks (Pelletier et al. 2019; Lemoine et al. 2019a,b), shear acceleration (Rieger & Duffy 2019; Tavecchio 2021), or magnetic reconnection (Blandford et al. 2017).

Many direct and indirect observations demonstrate the presence of a transverse stratified profile of AGN jets, characterized by the presence of a fast inner jet (spine) and a slower outer jet (sheath or layer), both inner and outer jets being relativistic. The most compelling observation of this structure at the limb-brightened jets was observed with very-long-baseline interferometry (VLBI), down to a few Schwarzschild radii for the nearest radio galaxies (Nagai et al. 2014; Kim et al. 2018). Seen at large angles, the slower layer is more Doppler boosted than the spine, leading to the appearance of a distinctive radio structure of an “empty pipe”. Observations of a different magnetic structure of the inner and outer jets via polarization measurements (Gabuzda et al. 2014; Avachat et al. 2016; Giroletti et al. 2004) support the idea that the two jet components may be launched from different processes, such as those proposed by Blandford & Znajek (1977) and Blandford & Payne (1982), for launching from the vicinity of the rotating supermassive black hole or from the accretion disk.

Theoretical studies in plasma physics also support this interpretation where the fast inner jet is responsible for most of the radiative output while having a lower density and a population dominated by electrons and positrons, while the outer jet is denser and less radiative, dominated by cold protons (Sol et al. 1989). The γ-ray detection of radio galaxies is challenging to explain when considering a uniform jet. It might be better explained by a stratified jet structure, where the particle and synchrotron fields of both jet components interact to produce a strong high-energy inverse-Compton emission (Tavecchio & Ghisellini 2008, 2014). The stratification of the relativistic jet can explain the spectral shape of the emission at multiwavelength, from the radio band to the X-ray band (Siemiginowska et al. 2007), and up to the (very) high energy γ-ray band (Ghisellini et al. 2005).

Large-scale magnetic fields seem to play an important role in extragalactic jets, in particular for the collimation of the jet and its acceleration (Porth et al. 2011; Fendt et al. 2012). Observations tend to show a correlation between the large-scale magnetic structure and the resulting synchrotron emission (Gabuzda et al. 2014; Walker et al. 2018).

The jet emission often shows a presence of bright spots (“knots”) that can be associated with standing and moving shocks. Such features have been detected in relativistic jets with radio and optical polarimetry observations (Perlman et al. 1999; Marshall et al. 2002), in the radio and millimeter band (Britzen et al. 2010; Jorstad et al. 2013) and in the X-ray band (Wilson & Yang 2002). One way to interpret these “knots” is by evoking recollimation shocks along the jet caused by pressure mismatches with the medium surrounding the jet (Marscher et al. 2008).

Flux variability is a characteristic feature of emission from radio-loud AGN and it depends on the observed frequency range and the AGN class. In the gigahertz radio bands, the shortest observed flare time scales can be of a few months in the case of BL Lac objects observed at a frequency of ν = 43 GHz (Wehrle et al. 2012, 2016), to a few years in the sample of several tens of radio-loud AGN observed at frequencies from ν = [4.8;  230] GHz (Hovatta et al. 2008; Nieppola et al. 2009). In the latter, the median flare duration was estimated to be 2.5 years, with flares occurring on average every four to six years. While in many cases the light curves of the detected flares exhibited a complex structure, sometimes including multiple peaks, in general, the decay times were found to be typically 1.3 times longer than the rise times. Their analysis concludes that the observed flare characteristics are in agreement with a generalized shock model (Valtaoja et al. 1992). In the case of the very nearby radio-galaxy M 87, VLBI observations (Walker et al. 2018) can locate fast flux variability from the radio core. At high energies, in X-rays and γ-rays, very rapid flares are observed from blazars and radio-galaxies with durations of days, down to time scales below an hour at teraelectronvolt energies.

Since the first analytical model (Blandford & Königl 1979) that was able to reproduce the flat radio spectra of jets with an isothermal jet associated with shock waves traveling in the jet, models have been evolving in several directions. Hydrodynamic and magneto-hydrodynamic (MHD) models are developed for in-depth studies of jet launching and propagation, while models of radiative transfer focus on the description of multiwavelength emission processes due to relativistic particle populations in the emission region, and particle in cell (PIC) models treat the microphysics of particle acceleration at small scales.

Today, increasingly sophisticated simulations address at the same time the macrophysics of the jet plasma and its radiative emission to try to improve our understanding of the jet physics from multiwavelength observations. Several hydrodynamical simulations of jets (Gómez et al. 1997; Agudo et al. 2001; Mimica et al. 2009) have shown that injection of perturbations at the base of the jet succeeds in reproducing the observed radio structure of superluminal and stationary components accounting for synchrotron radiation from a randomly oriented magnetic field. These simulations have also shown that perturbations traveling along an over-pressured jet can lead to the appearance of recollimation shocks.

Including a large-scale magnetic field structure in simulations of relativistic jets, Mizuno et al. (2015) studied the impact of the geometry of the magnetic field on recollimation shocks and rarefaction structures. They showed that the influence of the magnetic structure is not negligible and that, for example, axial fields lead to stronger collimation and rarefaction than toroidal fields. In studies by Martí (2015) and Martí et al. (2018), the authors simulated models of relativistic magnetized, axisymmetric jets with azimuthal velocity (i.e., rotation). For certain configurations, this azimuthal velocity leads to change the stationary shock-wave structure. Thus, they obtain a standing-shock structure and compute synthetic radio maps compatible with observations of parsec-scale extragalactic jets. Fuentes et al. (2018) are able to obtain polarization maps by computing the optically thin and linearly polarized synchrotron emission. They find that the electric vector polarization angles tend to be perpendicular to the axis of the jet between the recollimation shocks. This characteristic polarization can be compared with that obtained in VLBI observations of blazar jets.

Porth & Komissarov (2015) show that both unmagnetized and magnetized jets have great stability due to interactions with the ambient medium. The difference in pressure between the jet and the ambient medium allows the jet plasma to keep a causal connection. They propose an explanation for the Fanaroff-Riley dichotomy with different pressure values leading to the appearance of different structures of recollimation shocks. Simulations of stratified relativistic jets (Hervet et al. 2017) show that a two-component model of jets in interaction with an interstellar ambient medium can reproduce the observed knots through the generation of standing and moving shocks inside the jets.

Shock waves passing through a conical relativistic jet were first evoked by Marscher & Gear (1985) to interpret a flare of the quasar 3C 273 observed in 1983 from the millimeter to the infrared band. In this scenario, superluminal radio knots naturally arise as shocked regions in the jet. With today’s increased computing capacity, Fromm et al. (2016, 2018) have been able to reproduce typical flares observed at different wavelengths by simulating the interaction between ejecta and a recollimation shock structure.

Several characteristics of pc-scale relativistic jets are well reproduced with current models. Nevertheless, a better comprehension of the link between the supposed recollimation shock structure and the magnetic configuration is necessary to understand the multiwavelength observations, particularly during flares. We aim to understand the impact of the magnetic configuration in the jet on the dynamics of a perturbation (“ejecta”) at the base of the jet, which we suppose to be the cause for the observed flares. We have studied in detail its interaction with a two-component jet for different large-scale magnetic configurations, as well as the synchrotron emission during its propagation. The overall aim is to reproduce typical radio flare observation by the injection of such perturbations at the base of the jet and to put constraints on the magnetic field configuration and jet structure.

We carried out special relativistic magnetohydrodynamic (SRMHD) simulations of pc-scale jets with the MPI-AMRVAC code (Keppens et al. 2012), using a two-component model of relativistic jets with different magnetic configurations. Following Hervet et al. (2017) we considered an external jet component that carries more power than the internal one, while staying within the same order of magnitude. This kind of configuration leads to the formation of a typical standing-shock structure. We consider four different magnetic configurations (defined in cylindrical coordinates): hydrodynamic (H) with a turbulent magnetic field linked to the thermal pressure, toroidal (T) (with magnetic lines along φ), poloidal (P) (with magnetic lines along z, the jet propagation axis) and helical (HL) (with a magnetic pinch angle of 45°). Synchrotron radiation is computed in post-processing, assuming the injection of relativistic electrons following Gómez et al. (1995) and accounting for relativistic effects. In this way, radio light curves can be computed for the passage of ejecta in the stationary shock structure.

The organization of this paper is as follows. In Sect. 2, we briefly present the MPI-AMRVAC code and the numerical method used to solve the SRMHD equations, followed by a description of the numerical implementation of the two-component jet in Sect. 3. The structure of standing shocks that arises in the steady-state solution of this model is discussed in Sect. 4 for the four different magnetic-field configurations. The introduction of ejecta leads to the perturbations of the steady-state structure is developed in Sect. 5. The treatment of radiative transfer and generation of synchrotron maps and light curves in post-processing is explained in Sect. 6 and results are presented. To illustrate the relevance of our scenario to explain radio flares, recent results from observations of the blazar 3C 273 with the Very Long Baseline Array (VLBA) and Owens Valley Radio Observatory (OVRO) are discussed in Sect. 7 and are qualitatively compared to our simulations. Section 8 provides a general discussion of the implications of our scenario.

Throughout this paper we use natural units where the speed of light c = 1. As the distance unit will be the jet radius R_jet, the time unit will be R_jet in the co-mobile frame or R_jet/δ in the absolute frame (where δ is the Doppler factor).

2. Governing SRMHD equations and numerical method

We perform the numerical simulation of the relativistic magnetized two-component jet model using the 2D ideal special-relativistic magnetohydrodynamics version of the finite volume code MPI-AMRVAC in conservation form, using high-resolution shock-capturing methods (Meliani & Keppens 2007; Keppens et al. 2012). It solves the governing conservation equation as in Martí & Müller (2015) with U the state vector and Fⁱ the associated flux vectors:

$$\begin{array}{c}\hfill {\partial }_{\mathrm{t}}\mathit{U}+{\partial }_{x^i}{\mathit{F}}^i(\mathit{U})=0,\end{array}$$(1)

with:

$$\begin{array}{cc}& \mathit{U}={(D,S^j,\tau ,B^k)}^{\mathrm{T}},\hfill \end{array}$$(2)

$$\begin{array}{cc}& {\mathit{F}}^i=(D{v}^i,S^j{v}^i+p{\delta }^{\mathit{ij}}-b^j{B}^i/\gamma ,\tau v^i+p{v}^i-b^0{B}^i/\gamma ,v^i{B}^k\hfill \\ \hfill & -v^k{B}^i{)}^{\mathrm{T}},\hfill \end{array}$$(3)

where the rest-mass density is D, the momentum density in the j-direction S^j and the total energy density τ, calculated in the absolute frame. They are given by:

$$\begin{array}{cc}& D=\rho \gamma ,\hfill \end{array}$$(4)

$$\begin{array}{cc}& S^j=\rho h^{\ast }{\gamma }^2{v}^j-b^0{b}^j,\hfill \end{array}$$(5)

$$\begin{array}{cc}& \tau =\rho h^{\ast }{\gamma }^2-p_{\mathrm{tot}}-{(b^0)}^2,\hfill \end{array}$$(6)

where h^* = (1+ϵ_the+p_th/ρ+b²/ρ) with ϵ_the the internal energy. We note bⁱ the three vector magnetic field described in the co-moving frame (B and v describe the magnetic field and the velocity vector in the absolute frame) as:

$$\begin{array}{cc}& b^0=\gamma \mathit{B}.\mathit{v},\hfill \end{array}$$(7)

$$\begin{array}{cc}& b^i=B^i/\gamma +b^0{v}^i.\hfill \end{array}$$(8)

Finally, the Lorentz factor is given by $\gamma =1/\sqrt{1-{\mathit{v}}^i{\mathit{v}}_i}$ (with i running over the spatial indices [1, 2, 3]).

As a closure equation for the hydrodynamic part, we exploit the Synge equation of state (Mathews 1971; Meliani et al. 2004),

$$\begin{array}{c}\hfill p=\frac{(\mathrm{\Gamma }-1)\rho }{2}(\frac{ϵ}{\rho }-\frac{\rho }{ϵ}),\end{array}$$(9)

and the corresponding effective polytropic index is given as (Meliani et al. 2008),

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{eff}}=\mathrm{\Gamma }-{\displaystyle \frac{\mathrm{\Gamma }-1}{2}}(1-{\displaystyle \frac{1}{ϵ}}),\end{array}$$(10)

with ϵ = (1 + ϵ_the) is the specific internal energy. We fix Γ = 5/3; the effective index can vary in time and location between its relativistic (Γ = 4/3, when the thermal energy becomes larger than the mass energy) and its classical value (Γ = 5/3, when the thermal energy becomes negligible compared to the mass energy).

The divergence free condition for the magnetic field is satisfied by using the divergence cleaning method described by Dedner et al. (2002). We use the Harten–Lax–van Leer–Contact (HLLC) Riemann solver (Mignone & Bodo 2006) with third order reconstruction method cada3 (Cada & Torrilhon 2009). The combination of cada3 reconstruction (third order accurate) and HLLC flux computations is extremely robust and handles both sharp discontinuities and shock development accurately.

In the study of recollimation shocks, it is important to detect the shocks and distinguish them from compression waves in the numerical simulation. These internal shocks are in some cases very weak, making their detection difficult. For this purpose, in the hydrodynamics case, we use the shock detection algorithm by Zanotti et al. (2010). For the magnetized cases, we use a jump condition on the fast-magnetosonic number. We should note that in the simulation of the magnetized jet with only a toroidal-field component, the two methods converge at the vicinity of the jet axis.

3. Two-component model of a magnetized relativistic jet

In order to investigate the effect of the magnetic field configuration on the shock structure in a transverse structured jet, we use the two-component jet model proposed by Meliani & Keppens (2007, 2009) and Hervet et al. (2017). We adopt typical characteristics of a radio loud relativistic AGN jet, with a total kinetic luminosity of L_kin = 10⁴⁶ erg s⁻¹ (Ghisellini et al. 2014), and the jet radius taken to be R_jet = 0.1 pc at the parsec scale as observed in the jet of M 87 (Biretta et al. 2002). Concerning the internal jet structure, we assume an external jet radius R_out = R_jet and an internal jet radius R_in = R_out/3. We assume that the outer jet component carries an important fraction f_out = 0.75 of the total kinetic luminosity,

$$\begin{array}{c}\hfill L_{\mathrm{out},\mathrm{kin}}=f_{\mathrm{out}}{L}_{\mathrm{kin}}=({\gamma }_{\mathrm{out}}{h}_{\mathrm{out}}-1){\rho }_{\mathrm{out}}{\gamma }_{\mathrm{out}}\pi (R_{\mathrm{out}}^2-R_{\mathrm{in}}^2)v_{z,\mathrm{out}},\end{array}$$(11)

with the remaining L_in,  kin = (1−f_out) L_kin = 0.25 L_kin carried by the inner jet component.

For the simulations performed in this paper as initial condition, we set a non-rotating, superfast, magnetosonic cylindrical jet surrounded by a uniform, unmagnetized ambient medium with high rest-mass density. The rest-mass density and the total pressure ratio of the outer jet (ρ_out,p_out) to the ambient medium (ρ_am = 10³ cm⁻³, p_am = 1dyn cm⁻²) is (η_ρ=ρ_out/ρ_am = 10⁻², η_out,  p=p_{jet,  out}/p_am = 1.0). We assume a more over-pressured inner jet, with a lower rest-mass density and total pressure ratio of the inner jet (ρ_in, p_in) to the ambient medium (ρ_am, p_am) of (η_ρ=ρ_in/ρ_am = 10⁻³, η_in, p=p_in/p_am = 1.5) (Gómez et al. 1995, 1997). Moreover, the inner jet is assumed to be faster (γ_in = 10) than the outer jet (γ_out = 3) (Giroletti et al. 2004).

To investigate the effects of the magnetic field on the recollimation shocks and therefore on the evolution of the ejecta, we have considered different topologies: hydrodynamic (H) (as reference case), toroidal (T), axial (poloidal) (P) and helical (HL). The jet magnetization is set through the local maximum magnetization parameters at the inner and the outer jet component. The magnetization parameter is given by,

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{M}}=\frac{{\mathit{B}}^2+{(\mathit{v}·\mathit{B})}^2/2}{\rho h}·\end{array}$$(12)

In all magnetized cases, the magnetization parameter is set σ_M,  in = 10⁻³ for the inner jet component and σ_M,  out = 10⁻⁴ for the outer jet component, sufficiently low to allow the Fermi 1 acceleration mechanism to be efficient (Lemoine & Pelletier 2010; Sironi et al. 2013; Plotnikov et al. 2018). In all cases the relativistic jet is kinematically dominated.

In the poloidal field case (P), the magnetic field is uniform and parallel to the jet flow in each component, and (σ_M,  in,σ_M,  out) are constants. In the toroidal and helical cases (T, HL), we adopt the same profile as in Meliani & Keppens (2009),

$$\begin{array}{c}\hfill B_{\phi }=\{\begin{array}{ccc}& B_{\phi ,\mathrm{in},0}{\left(\frac{R}{R_{\mathrm{in}}}\right)}^{a_{\mathrm{in}}/2};\hfill & \hfill \mathrm{if}R<R_{\mathrm{in}},\\ \hfill & B_{\phi ,\mathrm{out},0}{\left(\frac{R}{R_{\mathrm{in}}}\right)}^{a_{\mathrm{out}}/2};\hfill & \hfill \mathrm{if}R\ge R_{\mathrm{in}},\end{array}\end{array}$$(13)

(B_{φ, in, 0}, B_{φ, out, 0}) are respectively the toroidal magnetic field strength of the inner and the outer jet component, at their interface and they are deduced from (Eq. (12)) and we fix the exponents (a_in,  a_out) = (0.5,   − 2) as in Meliani & Keppens (2009). It should be noted that the value of these exponents can have an influence on the resulting pattern of the recollimation shocks and the moving shock. Concerning the helical-field case (HL), we chose the same configuration as for the toroidal-field case (Eq. (13)) and a constant axial field strength with constant magnetic pitch angle θ_B = 45°.

Finally, the thermal pressure profile p_th(r) is deduced by assuming for each component a transverse equilibrium among the pressure gradient and Lorentz force following (Meliani & Keppens 2009),

$$\begin{array}{c}\hfill p_{\mathrm{th}}(r)=p_{\mathrm{tot}}-((1-v_z^2)·\left({\displaystyle \frac{1}{a_{\mathrm{in}}+0.5}}\right)·{\displaystyle \frac{B_{\phi }^2}{2}}+{\displaystyle \frac{B_z^2}{2}})·\end{array}$$(14)

The initial distributions of B_φ and p_mag in the different jet components can be seen in Fig. A.1. The smooth transition between the two components of the jet is imposed using a hyperbolic tangential function with an extension of R_in/2 on the density, and magnetization parameter σ_M. We note that all physical parameters are calculated with respect to the relativistic jet kinetic energy flux, radius, Lorentz factor and magnetization.

To make the simulation with high resolution and large spatial extend more tractable, we assume axisymmetry. The simulations are carried out in a domain size [R,  Z] = [16,  200] R_jet; we take a base resolution of 64 × 1000 and we allow for five levels of refinement, achieving an effective resolution of 1024 × 16 000. For the initial state, the jet inlet is set over the full jet axis Z.

4. Results from steady-state jet simulations

In the following, we present simulation results for the four magnetic configurations (H, T, P, and HL). The simulations are run until they reach a steady state with a stationary shock pattern appearing within the two components of the jet over the full simulation box (Fig. 1).

Fig. 1. Snapshot: different types of structured jets (H, T, P, and HL) without an injected perturbation. The propagation of the jet is going from left to right along the increasing value of Z. For each case, the density contour (in log-scale) is drawn on the bottom side and the thermal energy contour on the top side. Units on x and y-axis in Rjet unit.

4.1. Hydrodynamic case (H)

The over-pressured inner jet expands inducing the development of multiple steady conical recollimation shocks and rarefaction waves along the inner and the outer jet component (Fig. 1, top left). The high inertia ratio between the inner and outer jet component $({\gamma }_{\text{in}}^2{h}_{\text{in}}{\rho }_{\text{in}})/({\gamma }_{\text{out}}^2{h}_{\text{out}}{\rho }_{\text{out}})\simeq 42$ enhances the influence of the inner component on the outer component, even if the inner component carries only 25% of the total kinetic energy flux.

In the jet inlet at the inner/outer jet component interface, the lateral expansion of the over-pressured inner jet is associated with the development of conical rarefaction waves propagating in the inner jet component toward the jet axis, and conical shock waves propagating toward the jet edge in the outer jet. These propagating waves form an angle α = (1/ℳ) with the jet axis, where ℳ = γ v/(γ_sc_s) is the relativistic Mach number of the jet component in which the waves propagate, c_s is the sound speed and the associated Lorentz factor ${\gamma }_{\mathrm{s}}=1/\sqrt{1-c_{\mathrm{s}}^2}$ (for more details see Hervet et al. 2017). These waves produce a pattern of successive and near-equidistant standing recollimation shocks with separation distance of δZ_shock = 2 R_jetℳ.

In the inner jet component, the rarefaction wave propagates inward and forms an angle α_in ∼ 5.5° with the jet axis; when it reaches the jet axis, it is reflected outward under the same angle. At the interface between the inner and outer jet components, the wave is partially transmitted as a rarefaction wave in the outer jet with an angle α_out = arctan(1/ℳ_out) ≃ 7.5° and it is partially reflected as a shock wave toward the jet axis inside the inner jet. Each time the wave is partially reflected at the inner or the outer jet border, it changes the type from rarefaction to shock wave and vice versa. We can clearly see this structure in the evolution of the inner-jet Lorentz factor along the Z-direction (Fig. 2, blue curve). The partial transmission from the inner toward the outer jet dampens the wave and its intensity decreases with distance, whereas the wave intensity associated with the outer jet component increases with distance, since it accumulates the successive waves transmitted from the inner component (Fig. 1). Moreover, each time the waves from inner and outer jet interact, a partial wave reflection is produced toward the jet axis. A pattern of multiple waves arises with a wavelength depending on the Mach number of the inner and outer jet and on the jet radius.

Fig. 2. Profile of the mean Lorentz factor along the propagation axis Z and for a radius between R = [0.1,  0.23] Rjet (in the inner jet). We show the profiles of the 4 cases studied with (H) in blue, (T) in green, (P) in orange and (HL) in red without ejecta.

The transverse ram pressure produced by this reflected wave causes an expansion of the jet. The expansion of the inner jet component with its higher effective inertia is more pronounced. The radial expansion stabilizes fairly quickly and the jet opening angle becomes θ_jet ≃ 0.05° for the parameters chosen in our simulations (Fig. 3, blue line).

Fig. 3. Radius of the external jet (Rout) as a function of the distance along the jet axis. We represent, in different colors, the hydrodynamic (H, lowest curve), the toroidal (T), the poloidal (P) and the helical case (HL) of jet. The opening angle is deduced from the slope of a linear function fitted to these curves.

4.2. Toroidal case (T)

In the toroidal case, a large-scale azimuthal magnetic field is set up in the inner jet and in the outer jet with respectively B_in, φ = 50 mG and B_out, φ = 5 mG. Due to the over-pressured inner jet, a recollimation structure arises in the jet, with a “diamond” structure formed by compression and rarefaction waves in the two components (Fig. 1 (top right) or Fig. 4). Since the magnetic field strength chosen for this case is weak, the kinetic energy flux and the inertia ratio between the outer and the inner jet component remain of the same order as for the hydrodynamic case (H). However, there are a few differences between the toroidal and hydrodynamic cases.

Fig. 4. Zoom on the base of the toroidal jet (T). A standing shock structure appears on the pressure map. The R-axis and the Z-axis are given in Rjet unit. The white lines represent the rarefaction and compression wave shock fronts.

The high Lorentz factor associated with the inner jet component of at least γ = 10 induces a strong radial electric field $E_{\mathrm{r}}=-B_{\phi }·\sqrt{1-1/{\gamma }^2}$ that decreases the efficiency of the radial magnetic tension to collimate the jet. This efficiency decreases more in the rarefaction region where the Lorentz factor can reach a higher value of γ = 20 and the magnetic field strength decreases. As a result, the jet expands radially in these zones and the rarefaction zones become more elongated in the Z direction. The shock wave at the recollimation knots is dampened and appears closer to the jet axis compared to the hydrodynamic case. Therefore, the stationary shock wave decollimates the jet which expands with an opening angle two times higher than the hydrodynamic case with θ_jet ≃ 0.10° (Fig. 3, orange curve). Nevertheless, thanks to the magnetic tension, we recover a higher value of magnetic energy (∼B²) in the inner jet compared to an axial magnetic field. With distance from the jet inlet (Z > 100 R_jet), the strength of the shock wave and the associated rarefaction wave decreases. As a result, the ram-pressure applied by these waves weakens and the radial expansion of the jet slows down (Fig. 3, green curve). In this region, a radial instability grows in the jet inducing oscillations in the density and Lorentz factor.

4.3. Poloidal case (P)

Now we consider the case of an initial large-scale axial and uniform magnetic field in the inner-jet B_in,  p = 50 mG and in the outer-jet B_out,  p = 5 mG. As in the hydrodynamic case, the inner-jet component is over-pressured, which leads to multiple steady conical recollimation shocks and rarefactions waves that emerge at the jet inlet and propagate along the jet (Fig. 1, top left). At the jet inlet, the low magnetization induces only a weak difference with the hydrodynamic case. The distance between two recollimation shocks remains of the same order. The difference appears at a larger distance, where the magnetized jet starts to decollimate.

When the shock waves interact with the inner and the outer jet interface, the jet undergoes a weak transverse expansion. In this case, the jet expansion induced by the shock wave is stronger than in the hydrodynamics case, with a jet opening angle three-times larger θ_jet ≃ 0.17° (Fig. 3, orange curve). The successive shock waves push the magnetic field toward the axis, increasing the magnetic pressure. As a result, the jet decollimates. Moreover, the rarefaction regions exhibit a stronger radial expansion, inducing more efficient acceleration.

The interaction between the shock waves and the poloidal magnetic field induces radial instabilities that can be associated with the development of thin structures. By perturbing the jet, these instabilities develop themselves at the interface where the expansion is at a maximum and where the poloidal magnetic pressure in Z-direction increases. These instabilities become more pronounced with distance from the core and disturb the jet. As a result, the intensity of the shock and rarefaction waves decreases with distance in comparison with the hydrodynamic case. In addition, an expansion of the external jet toward the internal jet is well marked. This expansion, due to the heating of the external jet and the magnetic pressure along the Z-axis, tends to compress and recollimate the internal jet and thus modifies the structure of stationary shocks. It is this compression that tends to stretch the shock structure in the internal jet along the propagation axis.

4.4. Helical case (HL)

Finally, we consider a helical case with similar characteristics as the poloidal case, with a magnetic field strength of B = 50 mG in the inner jet and B = 5 mG for the outer jet. The only difference is the helical structure with a pitch angle between the azimuthal and axial magnetic direction fixed to 45°.

Overall, the recollimation shock structure is similar to the poloidal one (Fig. 1, bottom right). As for the poloidal case, the jet decollimation occurs after the second recollimation knot at a distance Z ≃ 50 R_jet from the jet inlet, and the jet opening angle tends to θ_jet ≃ 0.17° (Fig. 3, red curve). A small difference in the jet expansion appears between the poloidal and helical case at a large distance that results from the toroidal magnetic tension that tends to collimate the jet.

As in the poloidal case, radial instabilities develop due to the axial magnetic pressure increasing in the rarefaction region and they explain the strong variation of the Lorentz factor along the Z-direction in the inner jet (Fig. 2, red curve).

We observe again an expansion of the external jet toward the internal jet. The poloidal component of the magnetic field still involves compression of the internal jet and an elongation of the stationary shock structure in the direction of propagation. In addition and similar to the toroidal case, the toroidal component of the magnetic field implies a higher value of magnetic energy (compared to the case without toroidal component) in the inner jet.

5. Results from simulations of moving ejecta

A promising scenario to explain the observed flux variability in AGN jets is to consider the propagation of shock waves within the jet. In our models, the shock wave is caused by an initially over-dense (ρ_e = 10³ρ_in) spherical ejecta with radius R_e = R_jet/6 (half of the inner jet radius) and placed initially at the jet axis (R = 0) at the distance Z = R_e from the inner boundary. Its Lorentz factor is the same as the one of the inner jet γ_e = γ_in. With this configuration, the kinetic energy flux associated with the ejecta is 10⁴⁷ erg s⁻¹. We should note that the thermal energy of the ejecta is negligible in comparison to kinetic energy. All the time in this section is given in the co-mobile frame (δ = 1). The ejecta reaches the edge of the simulation box at time t = 200 R_jet, but the simulations run until t = 230 R_jet to cover the jet relaxation phase after the passage of the shock wave induced by the ejecta. The jet is presented with a moving shock wave in Fig. 5 corresponding respectively to the co-moving times 135 R_jet. Movies showing the evolution of the thermal energy and density during the propagation of the moving shock wave are available online.

Fig. 5. Snapshot: different types of structured jets (H, T, P, and HL) with an injected perturbation. The generated moving shock wave is located at ∼135 Rjet from the base. The propagation of the jet is going from left to right along the increasing value of Z. For each case, the density contour (in log-scale) is drawn on the bottom side and the thermal energy contour on the top side. Units on x and y-axis in Rjet unit.

5.1. Hydrodynamic case (H)

For the hydrodynamic case (H) (Fig. 5, top left), the ejecta remains well confined within the inner jet during a short period t < 5 R_jet until it starts interacting with the first standing shock. During this phase, the ejecta and thus the moving shock wave undergoes an adiabatic acceleration in the first rarefaction zone to reach a Lorentz factor γ_ms ≃ 14 before it collides with the first standing shock. As a result, the thermal energy of the ejecta increases and as it evolves in the rarefaction region, it is accelerated once more (Fig. 6, blue curve). In interactions between the ejecta and internal shocks, all the gained energy of the ejecta is transferred to shock waves that continue to propagate mainly within the inner jet with low transverse energy loss. Globally, the velocity of the ejecta will follow the profile of the Lorentz factor of the internal jet without seeing its velocity drop to the minimum value of γ = 11. As the ejecta is traveling through the jet, it will perturb momentarily the stationnary shock structure.

Fig. 6. Evolution of the Lorentz factor of the moving shock as a function of time in the co-moving frame. We represent, in different colors, the hydrodynamic (H), toroidal (T), poloidal (P) and the helical case (HL) of jet.

5.2. Toroidal case (T)

In the toroidal case (T) (Fig. 5, top right), the interaction of the moving shock wave with the standing shocks has some similarities with the previous case. In the first rarefaction zone, the ejecta accelerates adiabatically before it starts interacting with the first shock. The resulting moving shock wave is subject to a strong radial expansion after this first interaction and slows down close to its initial value of γ_ms = 10 at Z ∼ 50 R_jet. Then the moving shock sees its velocity increase between each recollimation shock (Fig. 6, green curve). Due to its initial higher inertia compared to the surrounding jet, the ejecta undergoes a stronger interaction when its crosses the recollimation zones. Therefore, the thermal pressure of the ejecta increases more than the ambient flow. Afterwards, in the rarefaction zone, the higher pressure and inertia of the shocked ejecta starts to behave as a fireball within the surrounding jet. As a result, the moving shock wave induced by the ejecta reaches a higher Lorentz factor than the surrounding jet.

As mentioned before (Sect. 4.2), the rarefaction zones are larger than in the other cases, especially after the fourth standing shock, where the shock wave is strongly accelerated to a value of γ_ms ≃ 30 at Z ∼ 125 R_jet. Then the moving shock interacts strongly with all the following standing shocks and causes them to oscillate along the jet axis (with a typical oscillation time close to ∼13 R_jet).

5.3. Poloidal case (P)

In the poloidal case (P) (Fig. 5, bottom left), the moving shock wave undergoes a strong acceleration to reach γ_ms ≃ 21 before it interacts with the first stationary features and slows down to γ_ms ≃ 13 at Z ∼ 24 R_jet (Fig. 6, orange curve.) In a second phase, the resulting moving shock propagates through the successive rarefaction zones where it accelerates and the compression zones where it decelerates. The mean Lorentz factor of the moving shock increases with distance until it reaches a stationary shock at the distance Z ∼ 125 R_jet. This acceleration results from the expansion of the inner jet component in this region. Beyond, the moving shock continues the propagation in the inner jet component which is compressed by the outer jet, and transverse instabilities grow along the stationary shock wave. In this region, the rarefaction zones are smaller and they are subject to multiple small scale stationary shocks. As a result, the moving shock wave decelerates to reach a Lorentz factor of 14. As in the toroidal case, beyond the fourth stationary shock, the passage of the moving shock induces an oscillation.

5.4. Helical case (HL)

In the helical case (HL) (Fig. 5, bottom right), the moving shock, after crossing the first internal shock wave at z = 5 R_jet, also undergoes a strong acceleration in the first rarefaction. Like in the poloidal case, this strong acceleration is followed by a strong interaction of the moving shock with stationary shocks, leading to the deceleration of the moving shock to (γ_ms = 10) (Fig. 6, red curve.) Beyond the fourth standing shock, the moving shock accelerates again to γ_ms ≃ 20.

6. Modeling the radiative transfer

6.1. Radiative processes

In a post-processing step using a separate code, we evaluate the synchrotron radiation of an electron population in the radio band and solve the radiative transfer equation along a given line of sight. We construct synchrotron emission maps to study the variation of the flux observed in different zones along the jet over time, as well as light curves of the spectral flux density integrated over the full simulated jet. In each cell, the relativistic electrons population is set with a power law, as expected for shock acceleration,

$$\begin{array}{c}\hfill N_{\mathrm{e}}(\gamma )\mathrm{d}\gamma =K{\gamma }^{-\mathrm{p}}\mathrm{d}\gamma ,\end{array}$$(15)

where γ_min < γ < γ_max. Following Gómez et al. (1995), we define the normalization coefficient as,

$$\begin{array}{c}\hfill K={\left[{\displaystyle \frac{e_{\mathrm{th},\mathrm{e}}(p-2)}{1-C_{\mathrm{E}}^{2-\mathrm{p}}}}\right]}^{\mathrm{p}-1}{\left[{\displaystyle \frac{1-C_{\mathrm{E}}^{1-\mathrm{p}}}{n_{\mathrm{e}}(p-1)}}\right]}^{\mathrm{p}-2},\end{array}$$(16)

where e_th, e = ϵ_ee_th is the fraction of thermal energy carried by the electrons, n_e = ϵ_e n is the fraction of the electron number density with ϵ_e = 0.1, p = 2.2 is the index of the power law and the coefficient C_E = γ_max/γ_min is set to 10³. As a simplification, we do not take into account radiative losses of the relativistic electrons and assume,

$$\begin{array}{c}\hfill {\gamma }_{\min }={\displaystyle \frac{e_{\mathrm{th},\mathrm{e}}}{n_{\mathrm{e}}}}{\displaystyle \frac{p-2}{p-1}}{\displaystyle \frac{1-C_{\mathrm{E}}^{1-\mathrm{p}}}{1-C_{\mathrm{E}}^{2-\mathrm{p}}}}·\end{array}$$(17)

In the present application, we are focusing only on the radio band, where the effect of radiative cooling is the smallest.

The specific intensity for each cell with index “i” is determined in the frame of a stationary observer at the location of the source (quantities in the co-moving frame are noted x′),

$$\begin{array}{c}\hfill I_{\nu ;\mathtt{i}}=I_{\nu ;\mathtt{i}-1}exp(-{\tau }_{\nu })+S_{\nu ;\mathtt{i}}(1-exp(-{\tau }_{\nu })),\end{array}$$(18)

where τ_ν is the optical depth due to synchrotron self-absorption and S_ν is the synchrotron source function.

Figure 7 shows schematically how the contribution of each cell along the line of sight to the total specific intensity I_ν is estimated. It should be noted that here we do not account for light crossing time effects, which can lead to a superposition of the signal from different emission regions along the jet due to the relativistic movement of the jet and the final speed of signal propagation. In the radio band, this effect is expected to be of minor importance (Chiaberge & Ghisellini 1999).

Fig. 7. Schematic representation of the resolution of the radiative transfer equation. We sum the different contributions along the line of sight.

The specific intensity I_ν (in [erg s⁻¹ cm⁻² Hz⁻¹ sterad⁻¹]) depends on the specific emission coefficient j_ν, the absorption coefficient α_ν and optical depth τ_ν. They are transformed to the observer frame at the location of the source for each cell,

$$\begin{array}{cc}& j_{\nu }={\delta }^2{j}_{\nu \prime },\hfill \end{array}$$(19)

$$\begin{array}{cc}& {\alpha }_{\nu }={\delta }^{-1}{\alpha }_{\nu \prime },\hfill \end{array}$$(20)

$$\begin{array}{cc}& {\tau }_{\nu }={\tau }_{\nu \prime }.\hfill \end{array}$$(21)

These transformations (Rybicki & Lightman 1979) depend on the Doppler factor δ = (γ(1 − βcos(θ_obs)))⁻¹ with $\gamma =1/\sqrt{1-{\mathit{v}}^2}$ the bulk Lorentz factor of the material in the cell and θ_obs is the angle between the direction of the jet axis and the line of sight. The different quantities (j_ν, α_ν, and τ_ν) are estimated in each cell following the approximations given by Katarzyński et al. (2001) which are appropriate for the cases studied here. The synchrotron flux in the observer frame on Earth is determined by,

$$\begin{array}{c}\hfill F_{\nu }=\frac{S_{\mathrm{e}}}{d_{\mathrm{l}}^2}(1+z)I_{\nu },\end{array}$$(22)

with d_l the luminosity distance (assuming H₀ = 70 km s⁻¹, Mpc⁻¹) and S_e the typical emission surface.

In our study, as an illustration, we choose the distance corresponding to M 87 (z = 0.00428). Finally, we obtain a 2D flux map of the synchrotron emission. To provide images that can be eventually compared with real VLBI images, we smooth the simulated images with a typical beamwidth obtained in the radio domain for M 87 close to 1.6 R_jet. To be able to distinguish between the emission from the jet and from the ejecta, we adjust an asymmetric 2D Gaussian distribution to the ejecta at each time step and extract the flux from the fitted (2σ) region (Fig. 9).

To add synchrotron emission in the hydrodynamical case (H), a passive turbulent magnetic field is added during the post-processing step. In this case, we assume that the magnetic energy density is a fraction ϵ_B of the thermal energy density,

$$\begin{array}{c}\hfill B_{\mathrm{turb}}=\sqrt{{ϵ}_B{ϵ}_{\mathrm{the}}},\end{array}$$(23)

with ϵ_B = 0.1.

6.2. Synthetic synchrotron images and light curves

Synthetic images are build by integrating the radiative transfer equation (Eq. (18)) along a line of sight, for all four cases (H, T, P, HL) and for one viewing angle θ_obs = 90° (Fig. 8). We show the light curves realized for the four cases tested at θ_obs = 90° (Fig. 9) and for the helical case at θ_obs = 2° , 15° (Fig. 10). The light curve is computed by integrating the flux on the full simulation box for the three observation angles. For an observation angle of θ_obs = 90°, it is also computed by integrating only the emission issued from the moving shock wave using the method described in Sect. 6.1.

Fig. 8. Snapshot: synchrotron emission map of the different types of jets (H, T, P, and HL) without ejecta and stationnary. Each map represent the flux intensity in mJy unit. The R and Z-axis are given in Rjet unit. These maps represent a resolution of the radiative transfer equation with an angle between the jet propagation axis and the line of sight equal to θobs = 90° and ν = 1 GHz.

The synthetic image of the stationary hydrodynamic jet (H) emission obtained for the observation angle θ_obs = 90° is presented in Fig. 8 (top). The emission is dominated by the stationary knots within the inner jet component. The intensity and the size of the emitting knots decrease with distance as a result of the slow damping of shock strength with distance (Sect. 3). The contribution from the external jet component to the emission is negligible, since the stationary shock wave within this jet component is weak.

The interaction between the moving shock wave, induced by the ejecta, with the stationary shocks rapidly increases the injection of accelerated particles and the associated synchrotron emission, causing a flare. The emission from the moving shock wave increases each time it passes through a standing knot. Moreover, after each collision, the stationary compression wave within the jet cannot ensure the radial cohesion of the shocked knot, thus the heated knot undergoes adiabatic expansion until it reaches the interface between the inner-outer jet. Afterwards the knot cools down and contracts. A remnant emission is associated to this adiabatic expansion. This illumination is visible on the jet light curve obtained with observation angle θ_obs = 90° (Fig. 9, top left). It contributes to the emission during the decay phase of the observed flares (Fig. 9, top left). Moreover, each time a knot is shocked by the moving shock wave, it starts to stretch and to oscillate along the jet axis. This behavior is associated with local changes of jet characteristics, such as the Mach number ℳ and the inner jet radius R_jet, in. As a result, the structure of the standing shock waves evolves over a long period. This variability is associated with an increase of the base emission of the jet compared to the stationary case.

Fig. 9. Light curve obtained by integrating the total synchrotron flux emitted from a simulation box with a size of [R = 8,  Z = 200] Rjet. The computation of the light curve is realized from the four different cases of jets (H, T, P and HL). The flux is integrated from t = 0 the time of the injection of the ejecta, until t ∼ 2300 Rjet/δ with δ(θobs = 90° ) = 0.1 in the absolute mobile frame and with ν = 1 GHz. We separate the total flux (orange) in two component: the jet (green) and the moving shock wave (blue). Movies showing the temporal evolution of the synchrotron flux map and associated light curves are available online.

The synthetic image of the stationary state for the toroidal case (T) is shown in the second panel of Fig. 8. Like in the hydrodynamic case (H), the emission comes mainly from the standing shocks in the inner jet, but the emission from the region between the shocks is stronger. Moreover, the stationary shocks appear more elongated along the jet and especially beyond the fourth knot. This behavior results from the toroidal magnetic field (Sect. 4.2).

The moving shock within the jet increases the emission, and the interaction with the stationary knots produces flares (Fig. 9, bottom left). In comparison with the hydrodynamic case (H), the strength of the flares is variable with time. The strongest flare occurs at a time t ∼ 1100 R_jet/δ when the moving shock crosses the fourth (and strongest) standing shock. This interaction leads to a strong deformation and oscillation of the knots along the jet axis. The following relaxation of the knots is associated with a slow flux decline in the light curve and produces an asymmetric flare. This asymmetry occurs in the strongest flares and is enhanced when a knot splits in two after interacting strongly with the moving shock.

The synthetic image for the poloidal case (P) is shown in the panel on Fig. 8, where the emission of the stationary shock structure is represented. The emission comes mainly from the part of the stationary shocks very near to the jet axis, where the shock is sufficiently strong. In comparison to cases H and T, the steady shock wave is weaker in the poloidal case. A complex emission structure is visible due to the magnetic pressure along the Z-axis and the expansion of the outer toward the inner jet (Sect. 4.3). In this case, the variations in flux observed between shock and rarefaction zones are weak. Due to the expansion of the outer jet, the emission structure is more elongated along the jet axis and therefore less pronounced.

As before in the P case, the moving shock induced by the ejecta triggers flares as it crosses the steady knots (Fig. 9, top right) but the efficiency of its interaction is weaker than in the H and T cases. Indeed, the poloidal magnetic field tends to damp the strength of the steady shock waves and it limits the transverse expansion of the moving shock wave. The resulting light curve is characterized by weak flares. Furthermore, when the moving shock interacts with a standing shock, instabilities form along the steady shock wave. These instabilities will propagate through the jet in the form of several moving shocks and thus locally increase slightly the synchrotron emission.

Finally, the synthetic images of the helical case (HL) can be see in the last panel of Fig. 8, where the emission of the stationary shock structure is represented. As in all the others cases, the emission comes mainly from the shock zones within the inner jet. However, the observed synchrotron flux is twice as high as in the other cases, due to the combination of the toroidal field effects with strong emission emanating from the center of the knots and the poloidal field effects with diffuse emission resulting from the compression of the inner jet by the outer component.

As in the previous cases, the moving shock will increase locally the synchrotron flux emitted and in particular when it crosses a stationary knot. But despite the fact that a toroidal magnetic field component is present, the shape of the light curve is very similar to the poloidal case (P) at θ_obs = 90° (Fig. 9, bottom right). As the emission is extended along the axis in the internal jet, variability is observed essentially from the inner jet itself due to the propagation of the moving shock. The presence of a poloidal field allows the development of transverse instabilities along steady shock waves. These instabilities will propagate in the inner jet in the form of moving shocks behind the principal moving shock and increase the emission from the jet. We notice that the impact of a magnetic tension due to the toroidal component is visible when the moving shock interacts with the last shock zones. In fact, the moving shock interacts more strongly with the last standing knot where the most pronounced flare occurs. As before, the interaction of the moving shock with standing knot induces an knot oscillation along the jet axis.

To evaluate the impact of a reduced viewing angle on the shape of the light curve, we show the results for the hydrodynamic configuration (H) (Fig. 10) for θ_obs = 15° and θ_obs = 2° for 1 GHz. Compared to the one obtained for θ_obs = 90°, they show the expected increase in the overall flux due to stronger Doppler beaming. The flares are shorter in duration due to the Doppler effect, but also due to self absorption. In certain cases, the dense ejecta partially hides the shocked knot behind it. For an observation angle θ_obs = 2°, the self absorption effect is very significant. The mean intensity decreases with distance, since, as the ejecta propagates within the jet, it hides a large number of stationary knots.

Fig. 10. Light curve obtain by integrating the total synchrotron flux emitted from a simulation box with a size of [R = 8,  Z = 200] Rjet. The computation of the light curve is realized for the hydrodynamic case (H) of jet. The flux is integrated from t = 0 the time of the injection of the ejecta, until respectively t ∼ (90,  13) Rjet/δ for δ(θobs = 15° ) ≃ 2.57 (top) and δ(θobs = 2° ) ≃ 18 (bottom) in the absolute frame and ν = 1 GHz. The gray dotted line represent the exit of ejecta from the simulation box.

7. Comparison with radio observations, the case of 3C 273

Our study shows the complex link between ejecta propagating in magnetized two-flow jets with stationary shocks, and its observed radio variability. All simulations in the present study consider a kinetic power of the outer jet larger than the one of the inner jet with an initial ratio of 3:1. The ratio of two-flow kinetic powers was proposed in Hervet et al. (2017) as a critical criterion discriminating between types of VLBI radio jets, which are themselves associated with spectral types of blazars (Hervet et al. 2016; Piner & Edwards 2018; Lister et al. 2019). Two-flow jets with kinetic powers within the same order of magnitude, such as the ones simulated in this study, were found to be the most similar to FSRQ-type blazars (flat spectrum radio quasar).

In order to compare the results of our simulations with an astrophysical case, we focus on the radio observations of one of the brightest and best monitored FSRQs over decades, 3C 273 (B1226+023). Its redshift of z = 0.1583 (Strauss et al. 1992) translates into a scaling factor of 2.73 pc mas⁻¹ considering H₀ = 70 km s⁻¹ Mpc⁻¹. 3C 273 displays a peculiar mix of fast moving and quasi-stationary radio knots. This hybrid radio kinematic behavior is most often observed in intermediate blazars (Low or Intermediate frequency-peaked BL Lacs) (Hervet et al. 2016). However 3C 273 significantly differs from these sources in the way that quasi-stationary knots are visible only during low-activity periods of the jet, not continuously.

For this study, we used 15.3 GHz observations from the VLBA, analyzed by the MOJAVE team up to August 2019 (Lister et al., in prep.). Most of these data up to December 2016 are publicly available from Lister et al. (2019). We combine this dataset with observations, at the same frequency, from the Californian 40 m single-dish OVRO telescope, which provides public light curves from a monitoring program of Fermi-LAT blazars (Richards et al. 2011). Our goal is to see, in a qualitative way, how the observed radio-VLBI ejecta influence the overall jet light curve observed with OVRO, as well as the luminosity evolution during their propagation.

We consider the period 2008−2019, overlapping both OVRO and MOJAVE observations with a dense VLBA monitoring. We specifically focus on four fast moving radio knots (k22, k31, k37, k36) and two quasi-stationary knots (k32, k35) which were observed during this period. All these observations are gathered in Fig. 11. The moving knot k39 is not considered in our study as it is only referenced by MOJAVE beyond the k32 zone, and we suspect several wrong identifications between k39, k32, and k35 from their referenced positions and luminosities.

Fig. 11. 3C 273 observed at 15.3 GHz. Top panel: distance to the core of radio knots analyzed by MOJAVE. We focus on the firmly identified components (in color). Straight lines are linear extrapolations of the moving knots based on their first 4 observations. Horizontal dashed lines show the mean position in the jet of the two observed quasi-stationary knots k32 and k35, with the gray band displaying the 1 sigma dispersion around the mean. Middle panel: radio jet light curve observed by OVRO. Bottom panel: measured flux of the radio core and moving knots. Dashed lines (extrapolation assuming smooth core emission) indicate that the observed core variability is actually due to the flux increase of the emerging moving knots when they are indiscernible from the core due to the limits of the VLBA angular resolution. Vertical lines show the most likely time when the moving knots pass through the stationary zone defined by k35 and the purple dashed line, with its associated uncertainty in gray.

In Fig. 11, we see that quasi-stationary knots appear during a specifically long quiescent state of the source of 2 years. This low activity period is marked by the long decrease of the radio flux seen by OVRO starting from 2011 up to 2013, and also corresponds to an absence of radio ejecta. This observation suggests that the radio-jet of 3C 273 presents at least two quasi-stationary knots (k32, k35) in its nominal (quiescent, non-perturbed) state.

In a few instances, the quasi-stationary knots k32 and k35 seem absent of the observations where one would expect to see them (k35 in two measurements between 2010 and 2012; k32 in one measurement before 2011). While these disappearances could be due to observing conditions or instrumental limitations, such as a jet locally outshined by the moving component; they could also highlight a relaxation time for the jet structure to return to its non-perturbed state after the passage of a moving knot. Linear extrapolations of the motions of the ejecta show that the jet radio luminosity starts to increase when ejecta are emerging from the radio core. For each new ejecta, the OVRO flux increase pursues up to the first standing knot marked by k35 (0.36 ± 0.08 mas to the core), and then it decreases while ejecta continue their propagation. VLBA observations of the flux from individual ejecta show a similar behavior as OVRO. The observed link between VLBI jet kinematics and the radio flux variability in 3C 273 leads to a consistent picture that is in agreement with the scenario we are proposing with our simulations.

Firstly, there is a systematic association between the passage of ejecta through the first standing knot k35 and a large flux increase of the overall jet radio luminosity. This crossing is the main phenomenon triggering the radio variability of the source. In addition, the flux increase of the ejecta up to k35 matches what is expected if k35 is a marker of a strong recollimation shock, as the ejecta should undergo a strong acceleration that enhances its Doppler factor by its passage through a large rarefaction zone before this recollimation shock.

Secondly, the ejecta enter in an uninterrupted cooling phase after their passage through k35. This is expected when considering that a rarefaction zone should follow the standing shock k35. The second stationary knot marked by k32 does not appear to play an important role for the variability. This suggests that this second recollimation shock is much less powerful than the first one.

Finally, the absence of flare can be linked with the passage of the ejecta through the radio core. This may suggest that the radio core marks the jet expanding funnel rather than a first stationary shock (Fig. 4a in Daly & Marscher 1988). However studying the variability at higher energy, outside the synchrotron-self-absorbed frequencies, would be necessary to confirm this assumption.

8. Discussion

8.1. Effects of jet stratification and magnetic field structure

As seen in previous studies, in hydrodynamic (Gómez et al. 1997; Fromm et al. 2018) and magnetized jets (Mizuno et al. 2015; Porth & Komissarov 2015; Fuentes et al. 2018), we find the well known diamond structure of standing shocks with a clear succession of shock and rarefaction zones in an over-pressured jet. A regular standing-shock structure with a linear intensity decrease over the distance from the base is observed. We found a combined effect of the large-scale magnetic configuration and jet stratification on the details of the shock structure. As observed in Hervet et al. (2017), the jet stratification induces a larger variety of standing-shock structures due to the interferences between the characteristics waves from the two layers of the jet. The toroidal magnetic field strengthens the standing shocks by inducing intense rarefaction regions, where the plasma is strongly accelerated (Fig. 1 (top right)). Stratification leads to the development of Rayleigh-Taylor instabilities at the outer-inner jet interface and along the standing shock wave as observed in Toma et al. (2017) in the hydrodynamic case. In the poloidal and helical cases, the magnetic pressure along the Z-axis amplifies these instabilities along the jet. They grow and heat both jet components. Within the inner jet, these instabilities interfere with the standing shocks and lead to a smoother structure and the appearance of a turbulent region at large distance (see Fig. 5 (bottom)).

In a transverse structured jet, the presence of a structured magnetic field amplifies the jet opening angle compared to the hydrodynamic case (Fig. 3). This is especially apparent in the poloidal and in the helical cases. Even if an outer jet layer shields the inner part from a homogeneous ambient medium (Porth & Komissarov 2015), the magnetic field modifies the topology of the characteristic waves of the fluid. The presence of a toroidal magnetic field component tends to limit this transverse expansion as observed in Mizuno et al. (2015) and at large distances in our simulations. On the other hand, the poloidal magnetic field component induces instabilities as we saw; these instabilities lead to a jet decollimation at medium and large distances. It was shown by Martí (2015) that introducing an azimuthal velocity component of the jet flow leads to a centrifugal force that further increases the jet opening. This effect is not currently treated in our models, but we expect that it may have an impact on the standing shock structure, which we will evaluate in a future study.

8.2. Ejecta and associated flares in the light curve

In all cases, the moving shock interacts with the successive standing rarefaction zones, where it is accelerated adiabatically, and collides with standing shocks, where it is heated and decelerated. These interactions lead to the appearance of flares in the light curves, due to thermal energy increase (cf. also Gómez et al. 1997; Fromm et al. 2016). The presence of a toroidal component of the magnetic field ensures the cohesion of the moving shock and of the standing shocks. The fact that the moving and standing shock zones are very compact is reflected in the emission of intense and clearly marked flares. We recover a similar behavior in the hydrodynamic case. On the contrary, where a poloidal field is present, the interaction between the moving shock wave and the more diffuse steady shock zones is weaker and its associated flares are less pronounced.

As we saw, the large variety of flares obtained is related to the intensity of the different knots, stemming from the combination of the characteristic waves of the plasma. The outer jet component allows interferences between the stationary shock waves in the inner jet and those of the outer jet. For certain conditions, the two stationary shock waves can combine and lead to a particularly strong emission. In the toroidal case, this effect leads to a strong standing shock, with an important rarefaction zone behind it (cf. the fourth standing shock in our simulations). This standing shock region is linked to the luminous flare emission in the light curve (Fig. 9). After this interaction, the ejecta is strongly accelerated in the rarefaction zone and will lose its cohesion due to fast adiabatic expansion.

8.3. Temporal flare structure

The simulated flares in H, T, and P cases are characterized by a temporal asymmetry with a fast rise and a slower decay phase, even though this is not always clearly visible due to the varying strength of the standing shocks and the limits in temporal resolution. When the moving shock interacts with the standing shocks, it heats and compresses them. Afterwards, the knots decompress. These interactions induce the formation of trailing recollimation shocks, already observed by Agudo et al. (2001) and Mimica et al. (2009), which will perturb standing knots along the jet axis and make them oscillate. This is associated with remnant emission from the shocked knots during their adiabatic cooling phase. This process is observed in all cases, but most clearly in the cases with pronounced interaction between the moving shock and the standing shock (cases H and T). Indeed, after the strong flare, we see in the toroidal case a delay between the emission of the ejecta and the emission of the jet (Fig. 9). This is the emission signature of the shocked knot, which causes a slight asymmetry. This additional radiation counterpart will tend to soften the slope after the ejecta has passed. This behavior is in accordance with the observed flare structure described by Hovatta et al. (2008) and Nieppola et al. (2009) and obtained numerically from over-pressured jets by Gómez et al. (1997) and Fromm et al. (2016).

8.4. Light curve with small observation angle

As the angle decreases to 15° or 2° (Fig. 10), the effect of the absorption of the moving shock along the line of sight becomes more important at low frequencies. Thus, the observed flare intensity, duration and asymmetry decrease as the occultation by the moving shock becomes more important. Moreover, after the flare associated with the interaction of the moving shock with the most intense steady knot, the flux intensity decreases. This is also observed by Fromm et al. (2016) but, in our case, we have not taken into account the effect of the “light travel delay” of photons emitted in different parts of the jet, which would lead to smoother light curves with a longer decay. This behavior is well distinguishable with viewing angle θ_obs = 2°.

8.5. Comparison with observation of 3C 273

With our shock propagation model, we arrive at a consistent, for now qualitative interpretation of the flaring behavior observed in 3C 273. Detailed modeling of the observed light curve is out of the scope of this work and will be addressed in a future publication. The number of recollimation shocks present in our simulated jets is much higher than the number of standing knots observed in the case of 3C 273, mostly due to the purely cylindrical shape of our jets. In real jets that are imperfectly aligned, the superposition of knots along the line of sight and radio opacity might also lead to an obscuration of standing knots for the observer. VLBI observations of the most nearby radio-galaxy (M 87, e.g., Asada et al. 2014) and of other AGN (Hervet et al. 2016) show the presence of multiple standing features that may be interpreted as series of recollimation shocks. Before the moving shock waves interact with the first standing shock, extrapolation of the emitted flux seems to indicate that they are undergoing a strong acceleration. This phenomenon appears clearly in our results, with the moving shock passing through the first rarefaction zone before interacting with the first knot. The acceleration seems especially pronounced in the magnetized cases as shown in the Fig. 6.

After the passage of moving shock waves, the disruption of the standing knots can be explained by the dynamics and/or by the radiative transfer within the jet. Concerning the dynamics aspect, the trailing recollimation shocks can perturb momentarily the standing-shock structure, as we saw in our simulations. On the other hand, the apparent “disappearance” of standing knots, if confirmed, may suggest that they are simply obscured by the brighter moving knots (k31 or k37) that hide the quasi-stationary knot, like we see in our synthetic light curve.

In our model, the light curve (Fig. 12) obtained at frequency 15.3 GHz (OVRO/MOJAVE frequency) and viewing angle 2° (e.g., Hervet et al. 2016) shows the interaction of the moving shock with the first two stationary shocks. The light curve integrates the flux emitted by the whole jet like a single dish telescope. We can thus compare our results with the OVRO data, in particular with the k37 event, which is isolated from other events.

Fig. 12. Light curve obtain by integrating the total synchrotron flux emitted during the first three interactions between the moving shock and the standing shocks. The computation of the light curve is realized for the four different cases of jets (H, T, P, and HL) for δ(θobs = 2° ) ≃ 18 and ν = 15.3 GHz. The red dots represent the estimated beginning, maximum and end of the flare.

We should note that the observations from Jorstad et al. (2017) found a different viewing angle of 6.4 ± 2.4°. To investigate the effects of the viewing angle, which is not precisely known, we also compute the light curves for 6° (Fig. 13).

Fig. 13. Light curve obtained as in Fig. 12. The computation of the light curve is realized for the cases where significant flares were present (P and HL) for δ(θobs = 6° ) ≃ 10 and ν = 15.3 GHz. The red dots represent the estimated beginning, maximum and end of the flare.

In the simulations, as in the observations, we find a flare during each interaction between the moving shock wave and a standing shock. However, the typical flare duration is very different (∼800 days for k37 and ∼150 days for our flares on average). This could be due to differences in the morphology of the jet, the size of the ejecta and stationary knots, and the uncertain value of the observation angle of 3C 273.

The shape of the observed flares seems to show a small asymmetry. To quantify this effect, we used the method proposed by Roy et al. (2018) and Nalewajko (2013) to compare the doubling (or halving) time in the rise (or fall) phase of the flare. Applying the method to the k37 flare, we found ξ_k37 = 0.12 ± 0.03 where the fall time is superior than the rise time. Applying the same method to the simulated flares, we found respectively ξ_H = 0.14 ± 0.01, ξ_T = 0.13 ± 0.02, ξ_P = 0.12 ± 0.03 and ξ_HL = −0.29 ± 0.06 for the hydrodynamic, toroidal, poloidal and helical case for θ_obs = 2° (Fig. 12).

At θ_obs = 6°, the shape of the light curve changes significantly due to beaming effects, such that the partial superposition of the first three flares does not allow us to clearly determine the presence or absence of asymmetries. The peak of the second flare is only well visible in the P and HL cases (Fig. 13). In the HL case, no asymmetry is found, while in the P case it has switched signs compared to the simulation at θ_obs = 2° (we found for this case ξ_P = −0.23 ± 0.03). We should also note that the current version of our model does not take into account time delay effects, which may play an important role at small angles. A more detailed study is beyond the scope of this work and will be applied in a future study to dedicated simulations for a given data set.

The amplitude of the flares is on average much larger in 3C 273 where the variability in radio can reach around twice the baseline value, compared to an increase of ∼15% in the same-frequency, small angle simulations (Fig. 12). This difference depends on the same characteristics as the flare duration. For the general picture, the main difference with 3C 273 is the presence of only two visible recollimation shocks, where only the first one is linked to strong flares, compared to a greater number of standing shocks and resulting flares in our simulations. The number of knots observed is however strongly linked to the angular resolution and sensitivity of the VLBA. Another stationary knot very close to the core at ∼0.10−0.16 mas, noted A1, ST1, or S1 has been detected when observing the jet at 43.2 GHz (Jorstad et al. 2005; Savolainen et al. 2006; Lisakov et al. 2017; Kim et al. 2020).

As discussed in Sect. 8.2 for the simulations, this effect may be explained in 3C 273 by a strong interaction between shock waves in the inner jet and outer jet occurring at the position of k35, at 0.36 mas to the core. This zone can lead to a major outburst when interacting with a moving shock. It also has the specificity of disrupting the downstream shock structure, which would explain the weak presence of the downstream stationary shock k32, as well as the absence of significant flare event associated with it, and the apparent adiabatic cooling of ejecta moving outside of k35.

To directly model observed flares with our radiative SRMHD code, a more important opening of the jet, reflecting changes in the density profile of the ambient medium, will be required. This should lead to a shock structure dominated by a few standing shocks close to the core. An implementation of the “light travel delay” effect (cf. Chiaberge & Ghisellini 1999) and of radiative cooling will also be needed for a more realistic description of the radiative emission.

9. Conclusions

We have investigated the effect of the large-scale magnetic field on the standing shocks and their interaction with ejecta within a two-component relativistic jet. The associated light curves, which were computed at two radio frequencies (ν = (1,  15.3) GHz) and for several observation angles (θ_obs = (2° , 15° , 90° )), show a variety of flares with varying durations and amplitudes.

Two-component magnetized jets are characterized by a complex standing-shock structure due to the interaction of characteristic waves propagating in the two jet components. In this way, jet stratification leads to the appearance of radio knots with a range of intensities along the jet. This is especially apparent in the toroidal case, where we recover a strong standing shock, giving rise to a pronounced flare in the interaction with the moving shock wave. Temporal asymmetry associated to the relaxation phase of the shocked standing knot is well visible for the strongest flares. The introduction of large-scale magnetic fields is seen to cause an intrinsic opening of the jet with an opening angle up to three times larger than the hydrodynamic case for our jet configuration.

Our scenario of moving ejecta interacting with standing shocks inside a two-component jet provides a good description of the kinematics and light curves seen in the jet of the FSRQ type blazar 3C 273 with VLBI and single-dish radio observations. In a preliminary study, at observation angle of 2°, an asymmetry in the simulated flare profiles, with a fall-time that is longer than the rise-time, was seen for the hydrodynamic, toroidal and poloidal cases, consistent with what is observed in the OVRO data for this source.

Movies

Movies associated with Fig. 5:

Movie 1 of Fig. 5 (helico_thermal_density) Access here

Movie 2 of Fig. 5 (hydro_thermal_density) Access here

Movie 3 of Fig. 5 (polo_thermal_density) Access here

Movie 4 of Fig. 5 (toro_thermal_density) Access here

Movies associated with Fig. 9:

Movie 5 of Fig. 9 (helico_map_CL) Access here

Movie 6 of Fig. 9 (hydro_map_CL) Access here

Movie 7 of Fig. 9 (polo_map_CL) Access here

Movie 8 of Fig. 9 (toro_map_CL) Access here

Acknowledgments

The authors thank the anonymous reviewer for his peer review. The computations of the SRMHD results were carried out on the OCCIGEN cluster at CINES (https://www.cines.fr/) in Montpellier (project named lut6216, allocation A0050406842) and from the MesoPSL Cluster at PSL University (http://www.mesopsl.fr/) in the Observatory of Paris. This research has made use of data from the MOJAVE database that is maintained by the MOJAVE team (Lister et al. 2018), and from the OVRO 40 m monitoring program which was supported in part by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G, and NSF grants AST-0808050 and AST-1109911, and private funding from Caltech and the MPIfR. The authors wish to thank M. Lister for providing preliminary MOJAVE data points for this study. We wish also thank M. Lemoine for providing insightful discussion during the study.

References

Appendix A: Additional figure

Fig. A.1. Left: initial representation of the toroidal magnetic field strength in space following Eq. (13) and assuming (Bφ, in, 0, Bφ, out, 0) ≡ 1. Right: initial representation of the magnetic pressure following Eq. (14).

All Figures

	Fig. 1. Snapshot: different types of structured jets (H, T, P, and HL) without an injected perturbation. The propagation of the jet is going from left to right along the increasing value of Z. For each case, the density contour (in log-scale) is drawn on the bottom side and the thermal energy contour on the top side. Units on x and y-axis in Rjet unit.
In the text

	Fig. 2. Profile of the mean Lorentz factor along the propagation axis Z and for a radius between R = [0.1,  0.23] Rjet (in the inner jet). We show the profiles of the 4 cases studied with (H) in blue, (T) in green, (P) in orange and (HL) in red without ejecta.
In the text

	Fig. 3. Radius of the external jet (Rout) as a function of the distance along the jet axis. We represent, in different colors, the hydrodynamic (H, lowest curve), the toroidal (T), the poloidal (P) and the helical case (HL) of jet. The opening angle is deduced from the slope of a linear function fitted to these curves.
In the text

	Fig. 4. Zoom on the base of the toroidal jet (T). A standing shock structure appears on the pressure map. The R-axis and the Z-axis are given in Rjet unit. The white lines represent the rarefaction and compression wave shock fronts.
In the text

Fig. 5. Snapshot: different types of structured jets (H, T, P, and HL) with an injected perturbation. The generated moving shock wave is located at ∼135 Rjet from the base. The propagation of the jet is going from left to right along the increasing value of Z. For each case, the density contour (in log-scale) is drawn on the bottom side and the thermal energy contour on the top side. Units on x and y-axis in Rjet unit.

In the text

	Fig. 6. Evolution of the Lorentz factor of the moving shock as a function of time in the co-moving frame. We represent, in different colors, the hydrodynamic (H), toroidal (T), poloidal (P) and the helical case (HL) of jet.
In the text

	Fig. 7. Schematic representation of the resolution of the radiative transfer equation. We sum the different contributions along the line of sight.
In the text

	Fig. 8. Snapshot: synchrotron emission map of the different types of jets (H, T, P, and HL) without ejecta and stationnary. Each map represent the flux intensity in mJy unit. The R and Z-axis are given in Rjet unit. These maps represent a resolution of the radiative transfer equation with an angle between the jet propagation axis and the line of sight equal to θobs = 90° and ν = 1 GHz.
In the text

Fig. 9. Light curve obtained by integrating the total synchrotron flux emitted from a simulation box with a size of [R = 8,  Z = 200] Rjet. The computation of the light curve is realized from the four different cases of jets (H, T, P and HL). The flux is integrated from t = 0 the time of the injection of the ejecta, until t ∼ 2300 Rjet/δ with δ(θobs = 90° ) = 0.1 in the absolute mobile frame and with ν = 1 GHz. We separate the total flux (orange) in two component: the jet (green) and the moving shock wave (blue). Movies showing the temporal evolution of the synchrotron flux map and associated light curves are available online.

In the text

Fig. 10. Light curve obtain by integrating the total synchrotron flux emitted from a simulation box with a size of [R = 8,  Z = 200] Rjet. The computation of the light curve is realized for the hydrodynamic case (H) of jet. The flux is integrated from t = 0 the time of the injection of the ejecta, until respectively t ∼ (90,  13) Rjet/δ for δ(θobs = 15° ) ≃ 2.57 (top) and δ(θobs = 2° ) ≃ 18 (bottom) in the absolute frame and ν = 1 GHz. The gray dotted line represent the exit of ejecta from the simulation box.

In the text

Fig. 11. 3C 273 observed at 15.3 GHz. Top panel: distance to the core of radio knots analyzed by MOJAVE. We focus on the firmly identified components (in color). Straight lines are linear extrapolations of the moving knots based on their first 4 observations. Horizontal dashed lines show the mean position in the jet of the two observed quasi-stationary knots k32 and k35, with the gray band displaying the 1 sigma dispersion around the mean. Middle panel: radio jet light curve observed by OVRO. Bottom panel: measured flux of the radio core and moving knots. Dashed lines (extrapolation assuming smooth core emission) indicate that the observed core variability is actually due to the flux increase of the emerging moving knots when they are indiscernible from the core due to the limits of the VLBA angular resolution. Vertical lines show the most likely time when the moving knots pass through the stationary zone defined by k35 and the purple dashed line, with its associated uncertainty in gray.

In the text

	Fig. 12. Light curve obtain by integrating the total synchrotron flux emitted during the first three interactions between the moving shock and the standing shocks. The computation of the light curve is realized for the four different cases of jets (H, T, P, and HL) for δ(θobs = 2° ) ≃ 18 and ν = 15.3 GHz. The red dots represent the estimated beginning, maximum and end of the flare.
In the text

	Fig. 13. Light curve obtained as in Fig. 12. The computation of the light curve is realized for the cases where significant flares were present (P and HL) for δ(θobs = 6° ) ≃ 10 and ν = 15.3 GHz. The red dots represent the estimated beginning, maximum and end of the flare.
In the text

	Fig. A.1. Left: initial representation of the toroidal magnetic field strength in space following Eq. (13) and assuming (Bφ, in, 0, Bφ, out, 0) ≡ 1. Right: initial representation of the magnetic pressure following Eq. (14).
In the text