11institutetext: LUTH, Observatoire de Paris, CNRS, PSL, Université de Paris;
5 Place Jules Janssen, 92190 Meudon, France
11email<EMAIL_ADDRESS>22institutetext: Santa Cruz
Institute for Particle Physics and Department of Physics, University of
California at Santa Cruz, Santa Cruz,CA 95064, USA
22email<EMAIL_ADDRESS>
# Flux variability from ejecta in structured relativistic jets with large-
scale magnetic fields

G. Fichet de Clairfontaine Z. Meliani A. Zech 111111 O. Hervet 22

###### Abstract

Context. Standing and moving shocks in relativistic astrophysical jets are
very promising sites for particle acceleration to large Lorentz factors and
for the emission from the radio up to the $\gamma$-ray band. They are thought
to be responsible for at least part of the observed variability in radio-loud
active galactic nuclei.

Aims. We aim to simulate the interactions of moving shock waves with standing
recollimation shocks in structured and magnetized relativistic jets and to
characterize the profiles of connected flares in the radio light curve.

Methods. Using the relativistic magneto-hydrodynamic (MHD) code MPI-AMRVAC and
a radiative transfer code in post-processing, we explore the influence of the
magnetic-field configuration and transverse stratification of an over-
pressured jet on its morphology, on the moving shock dynamics, and on the
emitted radio light curve. First, we investigate different large-scale
magnetic fields with their effects on the standing shocks and on the
stratified jet morphology. Secondly, we study the interaction of a moving
shock wave with the standing shocks. We calculated the synthetic synchrotron
maps and radio light curves and analyze the variability at two frequencies 1
and 15.3 GHz and for several observation angles. Finally, we compare the
characteristics of our simulated light curves with radio flares observed from
the blazar 3C 273 with the Owens Valley Radio Observatory (OVRO) and Very Long
Baseline Array (VLBA) in the MOJAVE survey between 2008 and 2019.

Results. We find that in a structured over-pressured relativistic jet, the
presence of the large-scale magnetic field structure changes the properties of
the standing shock waves and leads to an opening in the jet. The interaction
between waves from inner and outer jet components can produce strong standing
shocks. When crossing such standing shocks, moving shock waves accompanying
overdensities injected in the base of the jet cause very luminous radio
flares. The observation of the temporal structure of these flares under
different viewing angles probes the jet at different optical depths. At 1 GHz
and for small angles, the self-absorption caused by the moving shock wave
becomes more important and leads to a drop in the observed flux after it
interacts with the brightest standing knot. A weak asymmetry is seen in the
shape of the simulated flares, resulting from the remnant emission of the
shocked standing shocks. The characteristics of the simulated flares and the
correlation of peaks in the light curve with the crossing of moving and
standing shocks favor this scenario as an explanation of the observed radio
flares of 3C 273.

###### Key Words.:

magnetohydrodynamics (MHD) – ISM: jets and outflows – radiation mechanisms:
nonthermal – galaxies: active – quasars: individual (3C 273) – methods :
numerical

##  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. 2019b, a), shear acceleration
(Rieger & Duffy 2019; Tavecchio 2020), 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 gamma-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 gamma-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 $\nu=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
$\nu=\left[4.8;\leavevmode\nobreak\ 230\right]\leavevmode\nobreak\ \text{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 gamma-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 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 $\varphi$), poloidal (P) (with magnetic lines along $z$, the jet
propagation axis) and helical (HL) (with a magnetic pinch angle of
$45\degree$). 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 Section 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 Section 3. The structure of standing shocks that arises
in the steady-state solution of this model is discussed in Section 4 for the
four different magnetic-field configurations. The introduction of ejecta leads
to the perturbations of the steady-state structure is developed in Section 5.
The treatment of radiative transfer and generation of synchrotron maps and
light curves in post-processing is explained in Section 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 VLBA
and OVRO are discussed in Section 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_{\rm jet}$, the time unit will be
$R_{\rm jet}$ in the co-mobile frame or $R_{\rm jet}/\delta$ in the absolute
frame (where $\delta$ 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 $\@vec{U}$ the state vector and
$\@vec{F}^{\rm i}$ the associated flux vectors :

$\partial_{\rm t}\,\@vec{U}+\partial_{\rm x^{i}}\,\@vec{F}^{i}(\@vec{U})=0\,,$
(1)

with :

$\@vec{U}=\left(D,\,S^{\rm j},\,\tau,\,B^{\rm k}\right)^{\rm T}\,,$ (2)

$\displaystyle\@vec{F}^{\rm i}{}$ $\displaystyle=\Big{(}Dv^{\rm i},\,S^{\rm
j}v^{\rm i}+p\delta^{\rm ij}-b^{\rm j}B^{\rm i}/\gamma,\,\tau v^{\rm
i}+pv^{\rm i}-b^{\rm 0}B^{\rm i}/\gamma,\,v^{\rm i}B^{\rm k}$ (3)
$\displaystyle-v^{\rm k}B^{\rm i}\Big{)}^{\rm T}\,,$

where the rest-mass density is $D$, the momentum density in the j-direction
$S^{\rm j}$ and the total energy density $\tau$, calculated in the absolute
frame. They are given by :

$\displaystyle D=\rho\gamma\,,$ (4) $\displaystyle S^{\rm j}=\rho
h^{*}\gamma^{2}v^{\rm j}-b^{\rm 0}b^{\rm j}\,,$ (5) $\displaystyle\tau=\rho
h^{*}\gamma^{2}-p_{\rm tot}-\left(b^{\rm 0}\right)^{2}\,,$ (6)

where $h^{*}\,=\,\left(1+\epsilon_{\rm the}+p_{\rm th}/\rho+b^{2}/\rho\right)$
with $\epsilon_{\rm the}$ the internal energy. We note $b^{\rm i}$ the three
vector magnetic field described in the co-moving frame ($\@vec{B}$ and
$\@vec{v}$ describe the magnetic field and the velocity vector in the absolute
frame) as :

$\displaystyle b^{\rm 0}=\gamma\@vec{B}.\@vec{v}\,,$ (7) $\displaystyle b^{\rm
i}=B^{\rm i}/\gamma+b^{\rm 0}v^{\rm i}\,.$ (8)

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

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

$p=\frac{\left(\Gamma-1\right)\rho}{2}\left(\frac{\epsilon}{\rho}-\frac{\rho}{\epsilon}\right),$
(9)

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

$\Gamma_{\rm
eff}=\Gamma-\dfrac{\Gamma-1}{2}\left(1-\dfrac{1}{\epsilon}\right)\,,$ (10)

with $\epsilon=\left(1+\epsilon_{\rm the}\right)$ is the specific internal
energy. We fix $\Gamma=5/3$; the effective index can vary in time and location
between its relativistic ($\Gamma=4/3$, when the thermal energy becomes larger
than the mass energy) and its classical value ($\Gamma=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); Hervet et al. (2017). We
adopt typical characteristics of a radio loud relativistic AGN jet, with a
total kinetic luminosity of $L_{\rm kin}=10^{46}\leavevmode\nobreak\
\text{erg.s}^{-1}$ (Ghisellini et al. 2014), and the jet radius taken to be
$R_{\rm jet}\,=\,0.1\,{\rm 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_{\rm out}\,=\,R_{\rm jet}$ and an internal jet
radius $R_{\rm in}=R_{\rm out}/3$. We assume that the outer jet component
carries an important fraction $f_{\rm out}=0.75$ of the total kinetic
luminosity,

$L_{\rm out,\,kin}=f_{\rm out}\,L_{\rm kin}=\left(\gamma_{\rm out}h_{\rm
out}-1\right)\rho_{\rm out}\gamma_{\rm out}\pi\left(R^{2}_{\rm out}-R^{2}_{\rm
in}\right)v_{\rm z,out}\,,$ (11)

with the remaining $L_{\rm in,\,kin}=\left(1-f_{\rm out}\right)\,L_{\rm
kin}=0.25\,L_{\rm 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 $\left(\rho_{\rm out},p_{\rm
out}\right)$ to the ambient medium $\left(\rho_{\rm
am}=10^{3}\leavevmode\nobreak\ \text{cm}^{-3},p_{\rm am}=1\leavevmode\nobreak\
\text{dyn.cm}^{-2}\right)$ is $\left(\eta_{\rm\rho}=\rho_{\rm out}/\rho_{\rm
am}\,=\,10^{-2},\eta_{\rm out,p}=p_{\rm jet,\,out}/p_{\rm am}\,=\,1.0\right)$.
We assume a more over-pressured inner jet, with a lower rest-mass density and
total pressure ratio of the inner jet $\left(\rho_{\rm in},p_{\rm in}\right)$
to the ambient medium $\left(\rho_{\rm am},p_{\rm am}\right)$ of
$\left(\eta_{\rm\rho}=\rho_{\rm in}/\rho_{\rm am}\,=\,10^{-3},\eta_{\rm
in,p}=p_{\rm in}/p_{\rm am}\,=\,1.5\right)$ (Gómez et al. 1995, 1997).
Moreover, the inner jet is assumed to be faster ($\gamma_{\rm in}=10$) than
the outer jet ($\gamma_{\rm 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,

$\sigma_{\rm
M}=\frac{\@vec{B}^{2}+\left(\@vec{v}\cdot\@vec{B}\right)^{2}/2}{\rho\,h}\,.$
(12)

In all magnetized cases, the magnetization parameter is set $\sigma_{\rm
M,\,in}=10^{-3}$ for the inner jet component and $\sigma_{\rm
M,\,out}=10^{-4}$ 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 $\left(\sigma_{\rm M,\,in},\sigma_{\rm
M,\,out}\right)$ are constants. In the toroidal and helical cases (T, HL), we
adopt the same profile as in Meliani & Keppens (2009),

$B_{\varphi}=\left\\{\begin{aligned}
&B_{\rm\varphi,in,0}\,\left(\frac{R}{R_{\rm\,in}}\right)^{a_{\rm
in}/2}&&;\text{if}\ R<R_{\rm\,in}\,,\\\
&B_{\rm\varphi,out,0}\,\left(\frac{R}{R_{\rm\,in}}\right)^{a_{\rm
out}/2}&&;\text{if}\ R\geq R_{\rm\,in}\,,\end{aligned}\right.$ (13)

$\left(B_{\rm\varphi,in,0},\,B_{\rm\varphi,out,0}\right)$ 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_{\rm in},\,a_{\rm 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 $\theta_{B}=45^{\degree}$.

Finally, the thermal pressure profile $p_{\rm th}(r)$ is deduced by assuming
for each component a transverse equilibrium among the pressure gradient and
Lorentz force following (Meliani & Keppens 2009),

$p_{\rm th}(r)=p_{\rm tot}-\left(\left(1-v_{\rm
z}^{2}\right)\cdot\left(\dfrac{1}{a_{\rm
in}+0.5}\right)\cdot\dfrac{B_{\varphi}^{2}}{2}+\dfrac{B_{\rm
z}^{2}}{2}\right)\,.$ (14)

The initial distributions of $B_{\varphi}$ and $p_{\rm mag}$ in the different
jet components can be seen in Fig. 14.

The smooth transition between the two components of the jet is imposed using a
hyperbolic tangential function with an extension of $R_{\rm\,in}/2$ on the
density, and magnetization parameter $\sigma_{\rm 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 $\left[R,Z\right]=\left[16,200\right]\leavevmode\nobreak\ R_{\rm jet}$;
we take a base resolution of $64\times 1000$ and we allow for five levels of
refinement, achieving an effective resolution of $1024\times 16000$. 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, 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).

Figure 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 $R_{\rm jet}$ unit.

Figure 2: Profile of the mean Lorentz factor along the propagation axis $Z$
and for a radius between $R=\left[0.1,0.23\right]\leavevmode\nobreak\ R_{\rm
jet}$ (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.

###  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 $\left(\gamma^{2}_{\rm in}h_{\rm in}\rho_{\rm
in}\right)/\left(\gamma^{2}_{\rm out}h_{\rm out}\rho_{\rm out}\right)\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 $\alpha=\arctan\left(1/{\cal
M}\right)$ with the jet axis, where ${\cal M}=\gamma\,v/(\gamma_{\rm s}c_{\rm
s})$ is the relativistic Mach number of the jet component in which the waves
propagate, $c_{\rm s}$ is the sound speed and the associated Lorentz factor
$\gamma_{\rm s}=1/\sqrt{1-c_{\rm 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 $\delta Z_{\rm
shock}=2R_{\rm jet}{\cal M}$.
In the inner jet component, the rarefaction wave propagates inward and forms
an angle $\alpha_{\rm in}\sim 5.5\degree$ 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 $\alpha_{\rm
out}=\arctan\left(1/{\cal M_{\rm out}}\right)\simeq 7.5\degree$ 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.
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 $\theta_{\rm jet}\simeq 0.05\degree$
for the parameters chosen in our simulations (Fig. 3, blue line).

###  4.2 Toroidal case (T)

Figure 3: Radius of the external jet ($R_{\rm out}$) 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. Figure 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 $R_{\rm jet}$ unit. The white lines
represent the rarefaction and compression wave shock fronts.

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_{\rm in,\varphi}=50$ mG
and $B_{\rm out,\varphi}=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.

The high Lorentz factor associated with the inner jet component of at least
$\gamma=10$ induces a strong radial electric field
$E_{r}=-B_{\varphi}\cdot\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 $\gamma=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 $\theta_{\rm
jet}\simeq 0.10\degree$ (Fig. 3, orange curve). Nevertheless, thanks to the
magnetic tension, we recover a higher value of magnetic energy ($\sim B^{2}$)
in the inner jet compared to an axial magnetic field. With distance from the
jet inlet ($Z>100\leavevmode\nobreak\ R_{\rm 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_{\rm in,\leavevmode\nobreak\ p}\,=\,50\,{\rm mG}$
and in the outer-jet $B_{\rm out,\leavevmode\nobreak\ p}\,=\,5\,{\rm 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 $\theta_{\rm jet}\simeq 0.17\degree$
(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\degree$.
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\simeq
50\leavevmode\nobreak\ R_{\rm jet}$ from the jet inlet, and the jet opening
angle tends to $\theta_{\rm jet}\simeq 0.17\degree$ (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 ($\rho_{\rm e}=10^{3}\rho_{\rm
in}$) spherical ejecta with radius $R_{\rm e}=R_{\rm jet}/6$ (half of the
inner jet radius) and placed initially at the jet axis ($R=0$) at the distance
$Z=R_{\rm e}$ from the inner boundary. Its Lorentz factor is the same as the
one of the inner jet $\gamma_{\rm e}=\gamma_{\rm in}$. With this
configuration, the kinetic energy flux associated with the ejecta is
$10^{47}\leavevmode\nobreak\ \text{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 ($\delta=1$).
The ejecta reaches the edge of the simulation box at time
$t=200\leavevmode\nobreak\ R_{\rm jet}$, but the simulations run until
$t=230\leavevmode\nobreak\ R_{\rm 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\leavevmode\nobreak\ R_{\rm jet}$.

###  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\leavevmode\nobreak\
R_{\rm 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
$\gamma_{\rm ms}\simeq 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 $\gamma=11$. As
the ejecta is traveling through the jet, it will perturb momentarily the
stationnary shock structure.

###  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 $\gamma_{\rm ms}=10$ at $Z\sim
50\leavevmode\nobreak\ R_{\rm 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 (Section 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 $\gamma_{\rm ms}\simeq 30$ at
$Z\sim 125\leavevmode\nobreak\ R_{\rm 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 $\sim 13R_{\rm
jet}$).

Figure 5: Snapshot: different types of structured jets (H, T, P and HL) with
an injected perturbation. The generated moving shock wave is located at $\sim
135R_{\rm jet}$ 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 $R_{\rm jet}$ unit.

Figure 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.3 Poloidal case (P)

In the poloidal case (P) (Fig. 5, bottom left), the moving shock wave
undergoes a strong acceleration to reach $\gamma_{\rm ms}\simeq 21$ before it
interacts with the first stationary features and slows down to $\gamma_{\rm
ms}\simeq 13$ at $Z\sim 24\leavevmode\nobreak\ R_{\rm 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\sim
125\leavevmode\nobreak\ R_{\rm 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\leavevmode\nobreak\ R_{\rm
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 ($\gamma_{\rm ms}=10$) (Fig. 6, red
curve.) Beyond the fourth standing shock, the moving shock accelerates again
to $\gamma_{\rm ms}\simeq 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,

$N_{\rm
e}\left(\gamma\right)\text{d}\gamma=K\gamma^{-\text{p}}\leavevmode\nobreak\
\text{d}\gamma\,,$ (15)

where $\gamma_{\rm min}<\gamma<\gamma_{\rm max}$. Following (Gómez et al.
1995), we define the normalization coefficient as,

$K=\left[\dfrac{e_{\rm th,e}\left(\text{p}-2\right)}{1-C_{\rm
E}^{2-\text{p}}}\right]^{\leavevmode\nobreak\ \text{p}-1}\leavevmode\nobreak\
\left[\dfrac{1-C_{\rm E}^{1-\text{p}}}{n_{\rm
e}\left(\text{p}-1\right)}\right]^{\leavevmode\nobreak\ \text{p}-2}\,,$ (16)

where $e_{\rm th,e}=\epsilon_{\rm e}e_{\rm th}$ is the fraction of thermal
energy carried by the electrons, $n_{\rm e}=\epsilon_{\rm e}\,n$ is the
fraction of the electron number density with $\epsilon_{\rm e}=0.1$,
$\text{p}=2.2$ is the index of the power law and the coefficient $C_{\rm
E}=\gamma_{\rm max}/\gamma_{\rm min}$ is set to $10^{3}$. As a simplification,
we do not take into account radiative losses of the relativistic electrons and
assume,

$\gamma_{\rm min}=\dfrac{e_{\rm th,e}}{n_{\rm
e}}\dfrac{\text{p}-2}{\text{p}-1}\dfrac{1-C_{\rm E}^{1-\text{p}}}{1-C_{\rm
E}^{2-\text{p}}}\,.$ (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^{\prime}$),

$I_{\nu;\,\rm\texttt{i}}=I_{\nu;\,\rm\texttt{i}-1}\exp{\left(-\tau_{\nu}\right)}+S_{\nu;\,\rm\texttt{i}}\left(1-\exp{\left(-\tau_{\nu}\right)}\right)\,,$
(18)

where $\tau_{\nu}$ is the optical depth due to synchrotron self-absorption and
$S_{\nu}$ is the synchrotron source function.
Fig. 7 shows schematically how the contribution of each cell along the line of
sight to the total specific intensity $I_{\nu}$ 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).

Figure 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_{\nu}$ $\left(\leavevmode\nobreak\
\text{in}\leavevmode\nobreak\ \left[\text{erg}\leavevmode\nobreak\
\text{s}^{-1}\leavevmode\nobreak\ \text{cm}^{-2}\leavevmode\nobreak\
\text{Hz}^{-1}\leavevmode\nobreak\ \text{sterad}^{-1}\right]\right)$ depends
on the specific emission coefficient $j_{\nu}$, the absorption coefficient
$\alpha_{\nu}$ and optical depth $\tau_{\nu}$. They are transformed to the
observer frame at the location of the source for each cell,

$\displaystyle j_{\nu}$ $\displaystyle=\delta^{2}\leavevmode\nobreak\
j_{\nu^{\prime}}\,,$ (19) $\displaystyle\alpha_{\nu}$
$\displaystyle=\delta^{-1}\leavevmode\nobreak\ \alpha_{\nu^{\prime}}\,,$ (20)
$\displaystyle\tau_{\nu}$ $\displaystyle=\tau_{\nu^{\prime}}\,.$ (21)

These transformations (Rybicki & Lightman 1979) depend on the Doppler factor
$\delta=\left(\gamma\left(1-\beta\cos{\left(\theta_{\rm
obs}\right)}\right)\right)^{-1}$ with $\gamma=1/\sqrt{1-v^{2}}$ the bulk
Lorentz factor of the material in the cell and $\theta_{\rm obs}$ is the angle
between the direction of the jet axis and the line of sight. The different
quantities ($j_{\nu}$, $\alpha_{\nu}$ and $\tau_{\nu}$) 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,

$F_{\nu}=\frac{S_{\rm e}}{d_{l}^{2}}\left(1+z\right)I_{\nu}\,,$ (22)

with $d_{\rm l}$ the luminosity distance (assuming
$H_{0}=70\leavevmode\nobreak\ \text{km.s}^{-1}\text{Mpc}^{-1}$) and $S_{\rm
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\leavevmode\nobreak\ R_{\rm 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\sigma$) 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 $\epsilon_{\rm
B}$ of the thermal energy density,

$B_{\text{turb}}=\sqrt{\epsilon_{\rm B}\leavevmode\nobreak\ \epsilon_{\rm
the}}\,,$ (23)

with $\epsilon_{\rm 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 $\theta_{\text{obs}}={\leavevmode\nobreak\ 90\degree}$ (Fig. 8). We show
the light curves realized for the four cases tested at $\theta_{\rm
obs}={\leavevmode\nobreak\ 90\degree}$ (Fig. 9) and for the helical case at
$\theta_{\rm obs}={2\degree,\leavevmode\nobreak\ 15\degree}$ (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
$\theta_{\text{obs}}=90\degree$, it is also computed by integrating only the
emission issued from the moving shock wave using the method described in
Section 6.1.

Figure 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 $R_{\rm jet}$ 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
$\theta_{\text{obs}}=90\degree$ and $\nu=1$ GHz.

The synthetic image of the stationary hydrodynamic jet (H) emission obtained
for the observation angle $\theta_{\rm obs}=90\degree$ 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
(Section 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 $\theta_{\rm
obs}=90\degree$ (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 $\cal{M}$ and the inner jet
radius $R_{\rm 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.

Figure 9: Light curve obtained by integrating the total synchrotron flux
emitted from a simulation box with a size of
$\left[\text{R}=8,\leavevmode\nobreak\ \text{Z}=200\right]R_{\rm jet}$. 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 $\text{t}=0$ the time of
the injection of the ejecta, until $t\sim 2300\leavevmode\nobreak\ R_{\rm
jet}/\delta$ with $\delta(\theta_{\rm obs}=90\degree)=0.1$ in the absolute
mobile frame and with $\nu=1$ GHz. We separate the total flux (orange) in two
component : the jet (green) and the moving shock wave (blue).

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 (Section 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\sim 1100R_{\rm jet}/\delta$ 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 (Section 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
$\theta_{\rm obs}=90\degree$ (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 $\theta_{\rm obs}=15\degree$ and $\theta_{\rm obs}=2\degree$ for 1 GHz.
Compared to the one obtained for $\theta_{\rm obs}=90\degree$, 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 $\theta_{\rm obs}=2\degree$, 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.

Figure 10: Light curve obtain by integrating the total synchrotron flux
emitted from a simulation box with a size of
$\left[\text{R}=8,\leavevmode\nobreak\ \text{Z}=200\right]R_{\rm jet}$. The
computation of the light curve is realized for the hydrodynamic case (H) of
jet. The flux is integrated from $\text{t}=0$ the time of the injection of the
ejecta, until respectively $t\sim\left(90,\leavevmode\nobreak\
13\right)\leavevmode\nobreak\ R_{\rm jet}/\delta$ for $\delta(\theta_{\rm
obs}=15\degree)\simeq 2.57$ (top) and $\delta(\theta_{\rm obs}=2\degree)\simeq
18$ (bottom) in the absolute frame and $\nu=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_{0}=70$ km.s-1.Mpc-1. 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
40m 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 Figure 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.

Figure 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 Figure 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\pm 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.

Figure 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 $\delta(\theta_{\rm
obs}=2\degree)\simeq 18$ and $\nu=15.3\leavevmode\nobreak\ \rm{GHz}$. The red
dots represent the estimated beginning, maximum and end of the flare. Figure
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
$\delta(\theta_{\rm obs}=6\degree)\simeq 10$ and $\nu=15.3\leavevmode\nobreak\
\rm{GHz}$. The red dots represent the estimated beginning, maximum and end of
the flare.

##  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); Nieppola
et al. (2009) and obtained numerically from over-pressured jets by Gómez et
al. (1997); Fromm et al. (2016).

###  8.4 Light curve with small observation angle

As the angle decreases to $15\degree$ or $2\degree$ (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
$\theta_{\rm obs}=2\degree$.

###  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\degree$ (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.
We should note that the observations from Jorstad et al. (2017) found a
different viewing angle of $6.4\pm 2.4\degree$. To investigate the effects of
the viewing angle, which is not precisely known, we also compute the light
curves for $6\degree$ (Fig.13).

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 ($\sim 800$ days for k37 and $\sim
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); 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 $\xi_{\rm
k37}=0.12\pm 0.03$ where the fall time is superior than the rise time.
Applying the same method to the simulated flares, we found respectively
$\xi_{\rm H}=0.14\pm 0.01$, $\xi_{\rm T}=0.13\pm 0.02$, $\xi_{\rm P}=0.12\pm
0.03$ and $\xi_{\rm HL}=-0.29\pm 0.06$ for the hydrodynamic, toroidal,
poloidal and helical case for $\theta_{\rm obs}=2\degree$ (Fig. 12).

At $\theta_{\rm obs}=6\degree$, 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 $\theta_{\rm
obs}=2\degree$ (we found for this case $\xi_{\rm P}=-0.23\pm 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 $\sim 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 $\sim 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 Section 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 (Chiaberge & Ghisellini 1999, cf.) 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 ($\nu=\left(1,\leavevmode\nobreak\
15.3\right)\leavevmode\nobreak\ {\rm GHz}$) and for several observation angles
($\theta_{\rm obs}=\left(2\degree,\leavevmode\nobreak\
15\degree,\leavevmode\nobreak\ 90\degree\right)$), 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\degree$, 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.

###### Acknowledgements.

The authors thank the anonymous reviewer for his peer review. The computations
of the SRMHD results were carried out on the OCCIGEN cluster at
CINES111https://www.cines.fr/ in Montpellier (project named lut6216,
allocation A0050406842) and from the MesoPSL Cluster at PSL
University222http://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 40m 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

  * Agudo et al. (2001) Agudo, I., Gómez, J.-L., Martí, J.-M., et al. 2001, The Astrophysical Journal, 549, L183
  * Asada et al. (2014) Asada, K., Nakamura, M., Doi, A., Nagai, H., & Inoue, M. 2014, ApJ, 781, L2
  * Avachat et al. (2016) Avachat, S. S., Perlman, E. S., Adams, S. C., et al. 2016, The Astrophysical Journal, 832, 3
  * Biretta et al. (2002) Biretta, J. A., Junor, W., & Livio, M. 2002, New A Rev., 46, 239
  * Blandford et al. (2017) Blandford, R., Yuan, Y., Hoshino, M., & Sironi, L. 2017, Space Sci. Rev., 207, 291
  * Blandford & Koenigl (1979) Blandford, R. D. & Koenigl, A. 1979, Astrophys. Lett., 20, 15
  * Blandford & Königl (1979) Blandford, R. D. & Königl, A. 1979, The Astrophysical Journal, 232, 34
  * Blandford & Payne (1982) Blandford, R. D. & Payne, D. G. 1982, MNRAS, 199, 883
  * Blandford & Znajek (1977) Blandford, R. D. & Znajek, R. L. 1977, MNRAS, 179, 433
  * Britzen et al. (2010) Britzen, S., Kudryavtseva, N. A., Witzel, A., et al. 2010, A&A, 511, A57
  * Cada & Torrilhon (2009) Cada, M. & Torrilhon, M. 2009, J. Comput. Physics, 228, 4118
  * Chiaberge & Ghisellini (1999) Chiaberge, M. & Ghisellini, G. 1999, MNRAS, 306, 551
  * Daly & Marscher (1988) Daly, R. A. & Marscher, A. P. 1988, ApJ, 334, 539
  * Dedner et al. (2002) Dedner, A., Kemm, F., Kröner, D., et al. 2002, Journal of Computational Physics, 175, 645
  * Fendt et al. (2012) Fendt, C., Porth, O., & Vaidya, B. 2012, Journal of Physics: Conference Series, 372, 012011
  * Fromm et al. (2016) Fromm, C. M., Perucho, M., Mimica, P., & Ros, E. 2016, A&A, 588, A101
  * Fromm et al. (2018) Fromm, C. M., Perucho, M., Porth, O., et al. 2018, A&A, 609, A80
  * Fuentes et al. (2018) Fuentes, A., Gómez, J. L., Martí, J. M., & Perucho, M. 2018, The Astrophysical Journal, 860, 121
  * Gabuzda et al. (2014) Gabuzda, D. C., Reichstein, A. R., & O’Neill, E. L. 2014, MNRAS, 444, 172
  * Ghisellini et al. (2005) Ghisellini, G., Tavecchio, F., & Chiaberge, M. 2005, A&A, 432, 401
  * Ghisellini et al. (2014) Ghisellini, G., Tavecchio, F., Maraschi, L., Celotti, A., & Sbarrato, T. 2014, Nature, 515, 376
  * Giroletti et al. (2004) Giroletti, M., Giovannini, G., Feretti, L., et al. 2004, The Astrophysical Journal, 600, 127
  * Gómez et al. (1997) Gómez, J. L., Martí, J. M., Marscher, A. P., Ibáñez, J. M., & Alberdi, A. 1997, The Astrophysical Journall, 482, L33
  * Gómez et al. (1995) Gómez, J. L., Martí, J. M., Marscher, A. P., Ibáñez, J. M., & Marcaide, J. M. 1995, The Astrophysical Journal, 449
  * Hervet et al. (2016) Hervet, O., Boisson, C., & Sol, H. 2016, A&A, 592, A22
  * Hervet et al. (2017) Hervet, O., Meliani, Z., Zech, A., et al. 2017, å, 606, A103
  * Hovatta et al. (2008) Hovatta, T., Nieppola, E., Tornikoski, M., et al. 2008, å, 485, 51
  * Jorstad et al. (2005) Jorstad, S. G., Marscher, A. P., Lister, M. L., et al. 2005, AJ, 130, 1418
  * Jorstad et al. (2017) Jorstad, S. G., Marscher, A. P., Morozova, D. A., et al. 2017, ApJ, 846, 98
  * Jorstad et al. (2013) Jorstad, S. G., Marscher, A. P., Smith, P. S., et al. 2013, ApJ, 773, 147
  * Katarzyński et al. (2001) Katarzyński, K., Sol, H., & Kus, A. 2001, å, 367, 809
  * Keppens et al. (2012) Keppens, R., Meliani, Z., van Marle, A., et al. 2012, Journal of Computational Physics, 231, 718 , special Issue: Computational Plasma Physics
  * Kim et al. (2020) Kim, D.-W., Trippe, S., & Kravchenko, E. V. 2020, A&A, 636, A62
  * Kim et al. (2018) Kim, J. Y., Krichbaum, T. P., Lu, R. S., et al. 2018, A&A, 616, A188
  * Lemoine & Pelletier (2010) Lemoine, M. & Pelletier, G. 2010, MNRAS, 402, 321
  * Lemoine et al. (2019a) Lemoine, M., Pelletier, G., Vanthieghem, A., & Gremillet, L. 2019a, Phys. Rev. E, 100, 033210
  * Lemoine et al. (2019b) Lemoine, M., Vanthieghem, A., Pelletier, G., & Gremillet, L. 2019b, Phys. Rev. E, 100, 033209
  * Lisakov et al. (2017) Lisakov, M. M., Kovalev, Y. Y., Savolainen, T., Hovatta, T., & Kutkin, A. M. 2017, MNRAS, 468, 4478
  * Lister et al. (2018) Lister, M. L., Aller, M. F., Aller, H. D., et al. 2018, ApJS, 234, 12
  * Lister et al. (2019) Lister, M. L., Homan, D. C., Hovatta, T., et al. 2019, ApJ, 874, 43
  * Marscher & Gear (1985) Marscher, A. P. & Gear, W. K. 1985, The Astrophysical Journal, 298, 114
  * Marscher et al. (2008) Marscher, A. P., Jorstad, S. G., D’Arcangelo, F. D., et al. 2008, Nature, 452, 966
  * Marshall et al. (2002) Marshall, H. L., Miller, B. P., Davis, D. S., et al. 2002, ApJ, 564, 683
  * Martí & Müller (2015) Martí, J. M. & Müller, E. 2015, Living Reviews in Computational Astrophysics, 1, 3
  * Martí (2015) Martí, J.-M. 2015, MNRAS, 452, 3106
  * Martí et al. (2018) Martí, J. M., Perucho, M., Gómez, J. L., & Fuentes, A. 2018, International Journal of Modern Physics D, 27, 1844011
  * Mathews (1971) Mathews, W. G. 1971, ApJ, 165, 147
  * Meliani & Keppens (2007) Meliani, Z. & Keppens, R. 2007, å, 475, 785
  * Meliani & Keppens (2009) Meliani, Z. & Keppens, R. 2009, ApJ, 705, 1594
  * Meliani et al. (2008) Meliani, Z., Keppens, R., & Giacomazzo, B. 2008, A&A, 491, 321
  * Meliani et al. (2004) Meliani, Z., Sauty, C., Tsinganos, K., & Vlahakis, N. 2004, A&A, 425, 773
  * Mignone & Bodo (2006) Mignone, A. & Bodo, G. 2006, MNRAS, 368, 1040
  * Mimica et al. (2009) Mimica, P., Aloy, M. A., Agudo, I., et al. 2009, The Astrophysical Journal, 696, 1142
  * Mizuno et al. (2015) Mizuno, Y., Gómez, J. L., Nishikawa, K.-I., et al. 2015, The Astrophysical Journal, 809, 38
  * Nagai et al. (2014) Nagai, H., Haga, T., Giovannini, G., et al. 2014, ApJ, 785, 53
  * Nalewajko (2013) Nalewajko, K. 2013, Monthly Notices of the Royal Astronomical Society, 430, 1324
  * Nieppola et al. (2009) Nieppola, E., Hovatta, T., Tornikoski, M., et al. 2009, The Astronomical Journal, 137, 5022
  * Pelletier et al. (2019) Pelletier, G., Gremillet, L., Vanthieghem, A., & Lemoine, M. 2019, Phys. Rev. E, 100, 013205
  * Perlman et al. (1999) Perlman, E. S., Biretta, J. A., Zhou, F., Sparks, W. B., & Macchetto, F. D. 1999, ApJ, 117, 2185
  * Piner & Edwards (2018) Piner, B. G. & Edwards, P. G. 2018, ApJ, 853, 68
  * Plotnikov et al. (2018) Plotnikov, I., Grassi, A., & Grech, M. 2018, MNRAS, 477, 5238
  * Porth et al. (2011) Porth, O., Fendt, C., Meliani, Z., & Vaidya, B. 2011, The Astrophysical Journal, 737, 42
  * Porth & Komissarov (2015) Porth, O. & Komissarov, S. S. 2015, MNRAS, 452, 1089
  * Richards et al. (2011) Richards, J. L., Max-Moerbeck, W., Pavlidou, V., et al. 2011, ApJS, 194, 29
  * Rieger & Duffy (2019) Rieger, F. M. & Duffy, P. 2019, ApJ, 886, L26
  * Roy et al. (2018) Roy, N., Chatterjee, R., Joshi, M., & Ghosh, A. 2018, MNRAS, 482, 743
  * Rybicki & Lightman (1979) Rybicki, G. B. & Lightman, A. P. 1979, Radiative processes in astrophysics
  * Savolainen et al. (2006) Savolainen, T., Wiik, K., Valtaoja, E., & Tornikoski, M. 2006, A&A, 446, 71
  * Siemiginowska et al. (2007) Siemiginowska, A., Stawarz, Ł., Cheung, C. C., et al. 2007, The Astrophysical Journal, 657, 145
  * Sironi et al. (2013) Sironi, L., Spitkovsky, A., & Arons, J. 2013, The Astrophysical Journal, 771, 54
  * Sol et al. (1989) Sol, H., Pelletier, G., & Asseo, E. 1989, MNRAS, 237, 411
  * Strauss et al. (1992) Strauss, M. A., Huchra, J. P., Davis, M., et al. 1992, ApJS, 83, 29
  * Tavecchio (2020) Tavecchio, F. 2020, arXiv e-prints, arXiv:2011.03264
  * Tavecchio & Ghisellini (2008) Tavecchio, F. & Ghisellini, G. 2008, Monthly Notices of the Royal Astronomical Society: Letters, 385, L98
  * Tavecchio & Ghisellini (2014) Tavecchio, F. & Ghisellini, G. 2014, MNRAS, 443, 1224
  * Toma et al. (2017) Toma, K., Komissarov, S. S., & Porth, O. 2017, MNRAS, 472, 1253
  * Valtaoja et al. (1992) Valtaoja, E., Lahteenmaki, A., & Terasranta, H. 1992, Astronomy & Astrophysics, Suppl., 95, 73
  * Walker et al. (2018) Walker, R. C., Hardee, P. E., Davies, F. B., Ly, C., & Junor, W. 2018, The Astrophysical Journal, 855, 128
  * Wehrle et al. (2016) Wehrle, A. E., Grupe, D., Jorstad, S. G., et al. 2016, The Astrophysical Journal, 816, 53
  * Wehrle et al. (2012) Wehrle, A. E., Marscher, A. P., Jorstad, S. G., et al. 2012, ApJ, 758, 72
  * Wilson & Yang (2002) Wilson, A. S. & Yang, Y. 2002, ApJ, 568, 133
  * Zanotti et al. (2010) Zanotti, O., Rezzolla, L., Del Zanna, L., & Palenzuela, C. 2010, A&A, 523, A8

##  Appendix A Additional Figures

(a)

(b)

Figure 14: Left: Initial representation of the toroidal magnetic field
strength in space following Eq. 13 and assuming
$\left(B_{\rm\varphi,in,0},\,B_{\rm\varphi,out,0}\right)\equiv 1$. Right:
Initial representation of the magnetic pressure following Eq. 14.