Flux variability from ejectas in structured relativistic jets with large-scale magnetic fields
Gaëtan Fichet de Clairfontaine, Zakaria Meliani, Andreas Zech, Olivier Hervet
AAstronomy & Astrophysics manuscript no. 2021_01_flare © ESO 2021January 26, 2021
Flux variability from ejecta in structured relativistic jets withlarge-scale magnetic fields
G. Fichet de Clairfontaine Z. Meliani A. Zech and O. Hervet LUTH, Observatoire de Paris, CNRS, PSL, Université de Paris; 5 Place Jules Janssen, 92190 Meudon, Francee-mail: [email protected] Santa Cruz Institute for Particle Physics and Department of Physics, University of California at Santa Cruz, Santa Cruz,CA 95064,USAe-mail: [email protected]
January 26, 2021
ABSTRACT
Context.
Standing and moving shocks in relativistic astrophysical jets are very promising sites for particle acceleration to large Lorentzfactors and for the emission from the radio up to the γ -ray band. They are thought to be responsible for at least part of the observedvariability in radio-loud active galactic nuclei. Aims.
We aim to simulate the interactions of moving shock waves with standing recollimation shocks in structured and magnetizedrelativistic 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 di ff erent large-scale magnetic fields withtheir e ff ects on the standing shocks and on the stratified jet morphology. Secondly, we study the interaction of a moving shock wavewith the standing shocks. We calculated the synthetic synchrotron maps and radio light curves and analyze the variability at twofrequencies 1 and 15.3 GHz and for several observation angles. Finally, we compare the characteristics of our simulated light curveswith 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 changesthe properties of the standing shock waves and leads to an opening in the jet. The interaction between waves from inner and outer jetcomponents can produce strong standing shocks. When crossing such standing shocks, moving shock waves accompanying overden-sities injected in the base of the jet cause very luminous radio flares. The observation of the temporal structure of these flares underdi ff erent viewing angles probes the jet at di ff erent optical depths. At 1 GHz and for small angles, the self-absorption caused by themoving 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 standingshocks 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 nu-clei (AGN) has been detected from the radio band up to theteraelectronvolt range and shows frequent episodes of high fluxvariability in many sources. The emission is generally ascribedto a population of relativistic particles inside the jet that interactwith magnetic and photon fields to produce nonthermal radiationover a large wavelength range. In blazar-type AGN, the emissionfrom the jet plasma is amplified in the observer frame by rela-tivistic beaming e ff ects as the jet axis is aligned with the line ofsight, while the observed emission is weaker in radio galaxies,where the jet is seen under a larger angle. The multiwavelengthemission of AGN requires a mechanism able to reaccelerate par-ticles as they travel along the jet (Blandford & Koenigl 1979).The processes that are most frequently proposed are accelera-tion from internal shocks (Pelletier et al. 2019; Lemoine et al. 2019b,a), shear acceleration (Rieger & Du ff y 2019; Tavecchio2020), or magnetic reconnection (Blandford et al. 2017).Many direct and indirect observations demonstrate the pres-ence of a transverse stratified profile of AGN jets, characterizedby the presence of a fast inner jet (spine) and a slower outerjet (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 interfer-ometry (VLBI), down to a few Schwarzschild radii for the near-est radio galaxies (Nagai et al. 2014; Kim et al. 2018). Seen atlarge angles, the slower layer is more Doppler boosted than thespine, leading to the appearance of a distinctive radio structureof an "empty pipe." Observations of a di ff erent magnetic struc-ture 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 launchedfrom di ff erent processes, such as those proposed by Blandford& Znajek (1977) and Blandford & Payne (1982), for launching Article number, page 1 of 19 a r X i v : . [ a s t r o - ph . H E ] J a n & A proofs: manuscript no. 2021_01_flare from the vicinity of the rotating supermassive black hole or fromthe accretion disk.Theoretical studies in plasma physics also support this in-terpretation where the fast inner jet is responsible for most ofthe radiative output while having a lower density and a popula-tion dominated by electrons and positrons, while the outer jet isdenser and less radiative, dominated by cold protons (Sol et al.1989). The gamma-ray detection of radio galaxies is challeng-ing to explain when considering a uniform jet. It might be bet-ter explained by a stratified jet structure, where the particle andsynchrotron fields of both jet components interact to producea strong high-energy inverse-Compton emission (Tavecchio &Ghisellini 2008; Tavecchio & Ghisellini 2014). The stratificationof the relativistic jet can explain the spectral shape of the emis-sion at multiwavelength, from the radio band to the X-ray band(Siemiginowska et al. 2007), and up to the (very) high energygamma-ray band (Ghisellini et al. 2005).Large-scale magnetic fields seem to play an important role inextragalactic jets, in particular for the collimation of the jet andits acceleration (Porth et al. 2011; Fendt et al. 2012). Observa-tions tend to show a correlation between the large-scale magneticstructure 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 movingshocks. Such features have been detected in relativistic jets withradio and optical polarimetry observations (Perlman et al. 1999;Marshall et al. 2002), in the radio and millimeter band (Britzenet al. 2010; Jorstad et al. 2013) and in the X-ray band (Wilson& Yang 2002). One way to interpret these "knots" is by evok-ing recollimation shocks along the jet caused by pressure mis-matches with the medium surrounding the jet (Marscher et al.2008).Flux variability is a characteristic feature of emission fromradio-loud AGN and it depends on the observed frequency rangeand the AGN class. In the gigahertz radio bands, the shortest ob-served flare time scales can be of a few months in the case ofBL Lac objects observed at a frequency of ν =
43 GHz (Wehrleet al. 2012; Wehrle et al. 2016), to a few years in the sampleof 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.5years, with flares occurring on average every four to six years.While in many cases the light curves of the detected flares exhib-ited a complex structure, sometimes including multiple peaks, ingeneral, the decay times were found to be typically 1.3 timeslonger than the rise times. Their analysis concludes that the ob-served flare characteristics are in agreement with a generalizedshock model (Valtaoja et al. 1992). In the case of the very nearbyradio-galaxy M 87, VLBI observations (Walker et al. 2018) canlocate fast flux variability from the radio core. At high energies,in X-rays and gamma-rays, very rapid flares are observed fromblazars and radio-galaxies with durations of days, down to timescales 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 withan isothermal jet associated with shock waves traveling in thejet, models have been evolving in several directions. Hydrody-namic and MHD models are developed for in-depth studies ofjet launching and propagation, while models of radiative transferfocus on the description of multiwavelength emission processesdue to relativistic particle populations in the emission region,and particle in cell (PIC) models treat the microphysics of parti-cle acceleration at small scales. Today, increasingly sophisticated simulations address at thesame time the macrophysics of the jet plasma and its radia-tive emission to try to improve our understanding of the jetphysics from multiwavelength observations. Several hydrody-namical simulations of jets (Gómez et al. 1997; Agudo et al.2001; Mimica et al. 2009) have shown that injection of perturba-tions at the base of the jet succeeds in reproducing the observedradio structure of superluminal and stationary components ac-counting for synchrotron radiation from a randomly orientedmagnetic field. These simulations have also shown that pertur-bations traveling along an over-pressured jet can lead to the ap-pearance of recollimation shocks.Including a large-scale magnetic field structure in simula-tions of relativistic jets, Mizuno et al. (2015) studied the impactof the geometry of the magnetic field on recollimation shocksand rarefaction structures. They showed that the influence of themagnetic structure is not negligible and that, for example, axialfields lead to stronger collimation and rarefaction than toroidalfields. In studies by Martí (2015) and Martí et al. (2018), the au-thors simulated models of relativistic magnetized, axisymmetricjets with azimuthal velocity (i.e., rotation). For certain config-urations, this azimuthal velocity leads to change the stationaryshock-wave structure. Thus, they obtain a standing-shock struc-ture and compute synthetic radio maps compatible with observa-tions of parsec-scale extragalactic jets. Fuentes et al. (2018) areable to obtain polarization maps by computing the optically thinand linearly polarized synchrotron emission. They find that theelectric vector polarization angles tend to be perpendicular to theaxis of the jet between the recollimation shocks. This character-istic polarization can be compared with that obtained in VLBIobservations of blazar jets.Porth & Komissarov (2015) show that both unmagnetizedand magnetized jets have great stability due to interactions withthe ambient medium. The di ff erence in pressure between the jetand the ambient medium allows the jet plasma to keep a causalconnection. They propose an explanation for the Fanaro ff -Rileydichotomy with di ff erent pressure values leading to the appear-ance of di ff erent structures of recollimation shocks. Simulationsof stratified relativistic jets (Hervet et al. 2017) show that a two-component model of jets in interaction with an interstellar ambi-ent medium can reproduce the observed knots through the gen-eration of standing and moving shocks inside the jets.Shock waves passing through a conical relativistic jet werefirst evoked by Marscher & Gear (1985) to interpret a flare of thequasar 3C 273 observed in 1983 from the millimeter to the in-frared band. In this scenario, superluminal radio knots naturallyarise as shocked regions in the jet. With today’s increased com-puting capacity, Fromm et al. (2016, 2018) have been able to re-produce typical flares observed at di ff erent wavelengths by sim-ulating the interaction between ejecta and a recollimation shockstructure.Several characteristics of pc-scale relativistic jets are wellreproduced with current models. Nevertheless, a better compre-hension of the link between the supposed recollimation shockstructure and the magnetic configuration is necessary to un-derstand the multiwavelength observations, particularly duringflares. We aim to understand the impact of the magnetic config-uration in the jet on the dynamics of a perturbation ("ejecta")at the base of the jet, which we suppose to be the cause for theobserved flares. We have studied in detail its interaction with atwo-component jet for di ff erent large-scale magnetic configura-tions, as well as the synchrotron emission during its propagation.The overall aim is to reproduce typical radio flare observation by Article number, page 2 of 19. Fichet de Clairfontaine, Z. Meliani, A. Zech, O.Hervet: Flux variability from ejecta in structured relativistic jets. the injection of such perturbations at the base of the jet and to putconstraints 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 ofrelativistic jets with di ff erent magnetic configurations. Follow-ing Hervet et al. (2017) we considered an external jet componentthat carries more power than the internal one, while stayingwithin the same order of magnitude. This kind of configurationleads to the formation of a typical standing-shock structure.We consider four di ff erent magnetic configurations (defined incylindrical coordinates) : hydrodynamic (H) with a turbulentmagnetic field linked to the thermal pressure, toroidal (T) (withmagnetic lines along ϕ ), poloidal (P) (with magnetic lines along z , the jet propagation axis) and helical (HL) (with a magneticpinch angle of 45 ◦ ). Synchrotron radiation is computed inpost-processing, assuming the injection of relativistic electronsfollowing Gómez et al. (1995) and accounting for relativistice ff ects. In this way, radio light curves can be computed for thepassage of ejecta in the stationary shock structure.The organization of this paper is as follows. In Section 2, webriefly present the MPI-AMRVAC code and the numerical methodused to solve the SRMHD equations, followed by a descriptionof the numerical implementation of the two-component jet inSection 3. The structure of standing shocks that arises in thesteady-state solution of this model is discussed in Section 4 forthe four di ff erent magnetic-field configurations. The introductionof ejecta leads to the perturbations of the steady-state structureis developed in Section 5. The treatment of radiative transferand 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 theVLBA and OVRO are discussed in Section 7 and are qualita-tively compared to our simulations. Section 8 provides a generaldiscussion of the implications of our scenario.Throughout this paper we use natural units where the speedof 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 theabsolute frame (where δ is the Doppler factor).
2. Governing SRMHD equations and numericalmethod
We perform the numerical simulation of the relativistic mag-netized two-component jet model using the 2D ideal special-relativistic magnetohydrodynamics version of the finite volumecode
MPI-AMRVAC in conservation form, using high-resolutionshock-capturing methods (Meliani & Keppens 2007; Keppenset al. 2012). It solves the governing conservation equation as in(Martí & Müller 2015) with U the state vector and F i the asso-ciated flux vectors : ∂ t U + ∂ x i F i ( U ) = , (1)with : U = (cid:16) D , S j , τ, B k (cid:17) T , (2) F i = (cid:16) Dv i , S j v i + p δ ij − b j B i /γ, τ v i + pv i − b B i /γ, v i B k − v k B i (cid:17) T , (3) where the rest-mass density is D , the momentum density in thej-direction S j and the total energy density τ , calculated in theabsolute frame. They are given by : D = ργ , (4) S j = ρ h ∗ γ v j − b b j , (5) τ = ρ h ∗ γ − p tot − (cid:16) b (cid:17) , (6)where h ∗ = (cid:16) + (cid:15) the + p th /ρ + b /ρ (cid:17) with (cid:15) the the internal en-ergy. We note b i the three vector magnetic field described in theco-moving frame ( B and v describe the magnetic field and thevelocity vector in the absolute frame) as : b = γ B . v , (7) b i = B i /γ + b v i . (8)Finally, the Lorentz factor is given by γ = / (cid:112) − v i v i (with irunning over the spatial indices [1 , , p = ( Γ − ρ (cid:32) (cid:15)ρ − ρ(cid:15) (cid:33) , (9)and the corresponding e ff ective polytropic index is given as(Meliani et al. 2008), Γ e ff = Γ − Γ − (cid:32) − (cid:15) (cid:33) , (10)with (cid:15) = (1 + (cid:15) the ) is the specific internal energy. We fix Γ = / ff ective index can vary in time and location between itsrelativistic ( Γ = /
3, when the thermal energy becomes largerthan the mass energy) and its classical value (
Γ = /
3, whenthe thermal energy becomes negligible compared to the massenergy).The divergence free condition for the magnetic field is sat-isfied by using the divergence cleaning method described byDedner et al. (2002). We use the Harten–Lax–van Leer–Contact(HLLC) Riemann solver (Mignone & Bodo 2006) with third or-der reconstruction method cada3 (Cada & Torrilhon 2009). Thecombination of cada3 reconstruction (third order accurate) andHLLC flux computations is extremely robust and handles bothsharp discontinuities and shock development accurately.In the study of recollimation shocks, it is important to detect theshocks and distinguish them from compression waves in the nu-merical simulation. These internal shocks are in some cases veryweak, making their detection di ffi cult. For this purpose, in thehydrodynamics case, we use the shock detection algorithm byZanotti et al. (2010). For the magnetized cases, we use a jumpcondition on the fast-magnetosonic number. We should note thatin the simulation of the magnetized jet with only a toroidal-fieldcomponent, the two methods converge at the vicinity of the jetaxis. Article number, page 3 of 19 & A proofs: manuscript no. 2021_01_flare
3. Two-component model of a magnetizedrelativistic jet
In order to investigate the e ff ect of the magnetic field configura-tion on the shock structure in a transverse structured jet, we usethe two-component jet model proposed by Meliani & Keppens(2007, 2009); Hervet et al. (2017). We adopt typical characteris-tics of a radio loud relativistic AGN jet, with a total kinetic lumi-nosity of L kin = erg.s − (Ghisellini et al. 2014), and the jetradius taken to be R jet = . R out = R jet andan internal jet radius R in = R out /
3. We assume that the outer jetcomponent carries an important fraction f out = .
75 of the totalkinetic luminosity, L out , kin = f out L kin = ( γ out h out − ρ out γ out π (cid:16) R − R (cid:17) v z , out , (11)with the remaining L in , kin = (1 − f out ) L kin = . L kin carriedby the inner jet component.For the simulations performed in this paper as initialcondition, we set a non-rotating, superfast, magnetosoniccylindrical jet surrounded by a uniform, unmagnetized am-bient medium with high rest-mass density. The rest-massdensity and the total pressure ratio of the outer jet ( ρ out , p out )to the ambient medium (cid:16) ρ am = cm − , p am = − (cid:17) is (cid:16) η ρ = ρ out /ρ am = − , η out , p = p jet , out / p am = . (cid:17) . Weassume a more over-pressured inner jet, with a lowerrest-mass density and total pressure ratio of the in-ner jet ( ρ in , p in ) to the ambient medium ( ρ am , p am ) of (cid:16) η ρ = ρ in /ρ am = − , η in , p = p in / p am = . (cid:17) (Gómez et al.1995; Gómez et al. 1997). Moreover, the inner jet is assumed tobe faster ( γ in =
10) than the outer jet ( γ out =
3) (Giroletti et al.2004).To investigate the e ff ects of the magnetic field on the recol-limation shocks and therefore on the evolution of the ejecta, wehave considered di ff erent topologies: hydrodynamic ( H ) (as ref-erence case), toroidal ( T ), axial (poloidal) ( P ) and helical ( HL ).The jet magnetization is set through the local maximum magne-tization parameters at the inner and the outer jet component. Themagnetization parameter is given by, σ M = B + ( v · B ) / ρ h . (12)In all magnetized cases, the magnetization parameter is set σ M , in = − for the inner jet component and σ M , out = − for the outer jet component, su ffi ciently low to allow the Fermi1 acceleration mechanism to be e ffi cient (Lemoine & Pelletier2010; Sironi et al. 2013; Plotnikov et al. 2018). In all cases therelativistic jet is kinematically dominated.In the poloidal field case ( P ), the magnetic field is uni-form and parallel to the jet flow in each component, and (cid:0) σ M , in , σ M , out (cid:1) are constants. In the toroidal and helical cases ( T , HL ), we adopt the same profile as in Meliani & Keppens (2009), B ϕ = B ϕ, in , (cid:32) RR in (cid:33) a in / ; if R < R in , B ϕ, out , (cid:32) RR in (cid:33) a out / ; if R ≥ R in , (13) (cid:16) B ϕ, in , , B ϕ, out , (cid:17) are respectively the toroidal magnetic fieldstrength of the inner and the outer jet component, at their in-terface and they are deduced from (Eq. 12) and we fix the ex-ponents ( a in , a out ) = (0 . , −
2) as in (Meliani & Keppens 2009).It should be noted that the value of these exponents can have aninfluence on the resulting pattern of the recollimation shocks andthe moving shock.Concerning the helical-field case ( HL ), we chose the sameconfiguration as for the toroidal-field case (Eq. 13) and aconstant axial field strength with constant magnetic pitch angle θ B = ◦ .Finally, the thermal pressure profile p th ( r ) is deduced by as-suming for each component a transverse equilibrium among thepressure gradient and Lorentz force following (Meliani & Kep-pens 2009), p th ( r ) = p tot − (cid:16) − v (cid:17) · (cid:32) a in + . (cid:33) · B ϕ + B . (14)The initial distributions of B ϕ and p mag in the di ff erent jetcomponents can be seen in Fig. A.1.The smooth transition between the two components of the jetis imposed using a hyperbolic tangential function with an exten-sion of R in / σ M .We note that all physical parameters are calculated with respectto the relativistic jet kinetic energy flux, radius, Lorentz factorand magnetization.To make the simulation with high resolution and large spatialextend more tractable, we assume axisymmetry. The simulationsare carried out in a domain size [ R , Z ] = [16 , R jet ; we takea base resolution of 64 × ff ective resolution of 1024 × Z .
4. Results from steady-state jet simulations
In the following, we present simulation results for the four mag-netic configurations ( H , T , P , HL ). The simulations are run untilthey reach a steady state with a stationary shock pattern appear-ing within the two components of the jet over the full simulationbox (Fig. 1). H ) The over-pressured inner jet expands inducing the developmentof multiple steady conical recollimation shocks and rarefactionwaves along the inner and the outer jet component (Fig. 1,top left). The high inertia ratio between the inner and outerjet component (cid:16) γ h in ρ in (cid:17) / (cid:16) γ h out ρ out (cid:17) (cid:39)
42 enhances theinfluence of the inner component on the outer component, evenif the inner component carries only 25% of the total kineticenergy flux.
Article number, page 4 of 19. Fichet de Clairfontaine, Z. Meliani, A. Zech, O.Hervet: Flux variability from ejecta in structured relativistic jets. F i g . : S n a p s ho t: d i ff e r e n tt yp e s o f s t r u c t u r e d j e t s ( H , T , P a nd H L ) w it hou t a n i n j ec t e dp e r t u r b a ti on . T h e p r op a g a ti ono f t h e j e ti s go i ng fr o m l e f tt o r i gh t a l ong t h e i n c r ea s i ngv a l u e o f Z . F o r eac h ca s e , t h e d e n s it y c on t ou r( i n l og - s ca l e ) i s d r a w non t h e bo tt o m s i d ea nd t h e t h e r m a l e n e r gy c on t ou r on t h e t op s i d e . U n it s onx a ndy - a x i s i n R j e t un it . Article number, page 5 of 19 & A proofs: manuscript no. 2021_01_flare
Fig. 2: Profile of the mean Lorentz factor along the propagation axis Z and for a radius between R = [0 . , . R jet (in the innerjet). 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 jet inlet at the inner / outer jet component interface, thelateral expansion of the over-pressured inner jet is associatedwith the development of conical rarefaction waves propagatingin the inner jet component toward the jet axis, and conical shockwaves propagating toward the jet edge in the outer jet. Thesepropagating waves form an angle α = arctan (1 / M ) with the jetaxis, where M = γ v / ( γ s c s ) is the relativistic Mach number ofthe jet component in which the waves propagate, c s is the soundspeed and the associated Lorentz factor γ s = / (cid:112) − c (formore details see (Hervet et al. 2017)). These waves produce apattern of successive and near-equidistant standing recollima-tion shocks with separation distance of δ Z shock = R jet M .In the inner jet component, the rarefaction wave propagatesinward and forms an angle α in ∼ . ◦ with the jet axis; when itreaches the jet axis, it is reflected outward under the same angle.At the interface between the inner and outer jet components, thewave is partially transmitted as a rarefaction wave in the outerjet with an angle α out = arctan (1 / M out ) (cid:39) . ◦ and it is partiallyreflected as a shock wave toward the jet axis inside the inner jet.Each time the wave is partially reflected at the inner or the outerjet border, it changes the type from rarefaction to shock waveand vice versa. We can clearly see this structure in the evolutionof the inner-jet Lorentz factor along the Z -direction (Fig. 2,blue curve). The partial transmission from the inner toward theouter jet dampens the wave and its intensity decreases withdistance, whereas the wave intensity associated with the outerjet component increases with distance, since it accumulates thesuccessive waves transmitted from the inner component (Fig.1). Moreover, each time the waves from inner and outer jetinteract, a partial wave reflection is produced toward the jet axis.A pattern of multiple waves arises with a wavelength dependingon the Mach number of the inner and outer jet and on the jetradius.The transverse ram pressure produced by this reflected wave causes an expansion of the jet. The expansion of the inner jetcomponent with its higher e ff ective inertia is more pronounced.The radial expansion stabilizes fairly quickly and the jet openingangle becomes θ jet (cid:39) . ◦ for the parameters chosen in oursimulations (Fig. 3, blue line). T ) Fig. 3: Radius of the external jet ( R out ) as a function of the dis-tance along the jet axis. We represent, in di ff erent colors, the hy-drodynamic ( H , lowest curve), the toroidal ( T ), the poloidal ( P )and the helical case ( HL ) of jet. The opening angle is deducedfrom the slope of a linear function fitted to these curves. Article number, page 6 of 19. Fichet de Clairfontaine, Z. Meliani, A. Zech, O.Hervet: Flux variability from ejecta in structured relativistic jets.
Fig. 4: Zoom on the base of the toroidal jet ( T ). A standing shockstructure appears on the pressure map.The R -axis and the Z -axisare given in R jet unit. The white lines represent the rarefactionand compression wave shock fronts.In the toroidal case, a large-scale azimuthal magnetic fieldis set up in the inner jet and in the outer jet with respectively B in ,ϕ =
50 mG and B out ,ϕ = H ). However, there are a few di ff erences betweenthe toroidal and hydrodynamic cases.The high Lorentz factor associated with the inner jet com-ponent of at least γ =
10 induces a strong radial electric field E r = − B ϕ · (cid:112) − /γ that decreases the e ffi ciency of the radialmagnetic tension to collimate the jet. This e ffi ciency decreasesmore in the rarefaction region where the Lorentz factor can reacha higher value of γ =
20 and the magnetic field strength de-creases. As a result, the jet expands radially in these zones andthe rarefaction zones become more elongated in the Z direction. The shock wave at the recollimation knots is dampened and ap-pears closer to the jet axis compared to the hydrodynamic case.Therefore, the stationary shock wave decollimates the jet whichexpands with an opening angle two times higher than the hydro-dynamic case with θ jet (cid:39) . ◦ (Fig. 3, orange curve). Neverthe-less, thanks to the magnetic tension, we recover a higher valueof magnetic energy ( ∼ B ) in the inner jet compared to an axialmagnetic field. With distance from the jet inlet ( Z > R jet ), thestrength of the shock wave and the associated rarefaction wavedecreases. As a result, the ram-pressure applied by these wavesweakens and the radial expansion of the jet slows down (Fig. 3,green curve). In this region, a radial instability grows in the jetinducing oscillations in the density and Lorentz factor. P ) Now we consider the case of an initial large-scale axial anduniform magnetic field in the inner-jet B in , p =
50 mG and inthe outer-jet B out , p = ff erence with the hydrodynamic case. The distance betweentwo recollimation shocks remains of the same order. The di ff er-ence appears at a larger distance, where the magnetized jet startsto decollimate.When the shock waves interact with the inner and the outer jetinterface, the jet undergoes a weak transverse expansion. In thiscase, the jet expansion induced by the shock wave is strongerthan in the hydrodynamics case, with a jet opening angle three-times larger θ jet (cid:39) . ◦ (Fig. 3, orange curve). The successiveshock waves push the magnetic field toward the axis, increasingthe magnetic pressure. As a result, the jet decollimates. More-over, the rarefaction regions exhibit a stronger radial expansion,inducing more e ffi cient acceleration.The interaction between the shock waves and the poloidal mag-netic field induces radial instabilities that can be associated withthe development of thin structures. By perturbing the jet, theseinstabilities develop themselves at the interface where the expan-sion is at a maximum and where the poloidal magnetic pressurein Z -direction increases. These instabilities become more pro-nounced with distance from the core and disturb the jet. As aresult, the intensity of the shock and rarefaction waves decreaseswith distance in comparison with the hydrodynamic case. In ad-dition, an expansion of the external jet toward the internal jet iswell marked. This expansion, due to the heating of the externaljet and the magnetic pressure along the Z -axis, tends to compressand recollimate the internal jet and thus modifies the structure ofstationary shocks. It is this compression that tends to stretch theshock structure in the internal jet along the propagation axis. HL ) Finally, we consider a helical case with similar characteristics asthe poloidal case, with a magnetic field strength of B =
50 mG inthe inner jet and B = ff erenceis the helical structure with a pitch angle between the azimuthaland axial magnetic direction fixed to 45 ◦ .Overall, the recollimation shock structure is similar to thepoloidal one (Fig. 1, bottom right). As for the poloidal case, thejet decollimation occurs after the second recollimation knot at adistance Z (cid:39) R jet from the jet inlet, and the jet opening angle Article number, page 7 of 19 & A proofs: manuscript no. 2021_01_flare tends to θ jet (cid:39) . ◦ (Fig. 3, red curve). A small di ff erence in thejet expansion appears between the poloidal and helical case at alarge distance that results from the toroidal magnetic tension thattends to collimate the jet.As in the poloidal case, radial instabilities develop due to theaxial magnetic pressure increasing in the rarefaction region andthey 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 theinternal jet. The poloidal component of the magnetic field stillinvolves compression of the internal jet and an elongation of thestationary shock structure in the direction of propagation. In ad-dition and similar to the toroidal case, the toroidal componentof the magnetic field implies a higher value of magnetic energy(compared to the case without toroidal component) in the innerjet.
5. Results from simulations of moving ejecta
A promising scenario to explain the observed flux variability inAGN jets is to consider the propagation of shock waves withinthe jet. In our models, the shock wave is caused by an initiallyover-dense ( ρ e = ρ in ) spherical ejecta with radius R e = R jet / R =
0) at the distance Z = R e from the inner boundary. ItsLorentz factor is the same as the one of the inner jet γ e = γ in .With this configuration, the kinetic energy flux associated withthe ejecta is 10 erg / s. We should note that the thermal energyof the ejecta is negligible in comparison to kinetic energy. Allthe time in this section is given in the co-mobile frame ( δ = t = R jet , but the simulations run until t = R jet to cover the jetrelaxation phase after the passage of the shock wave induced bythe ejecta. The jet is presented with a moving shock wave in Fig.5 corresponding respectively to the co-moving times 135 R jet . H ) For the hydrodynamic case ( H ) (Fig. 5, top left), the ejecta re-mains well confined within the inner jet during a short period t < R jet until it starts interacting with the first standing shock.During this phase, the ejecta and thus the moving shock waveundergoes an adiabatic acceleration in the first rarefaction zoneto reach a Lorentz factor γ ms (cid:39)
14 before it collides with thefirst standing shock. As a result, the thermal energy of the ejectaincreases and as it evolves in the rarefaction region, it is accel-erated once more (Fig. 6, blue curve). In interactions betweenthe ejecta and internal shocks, all the gained energy of the ejectais transferred to shock waves that continue to propagate mainlywithin the inner jet with low transverse energy loss. Globally,the velocity of the ejecta will follow the profile of the Lorentzfactor of the internal jet without seeing its velocity drop to theminimum value of γ =
11. As the ejecta is traveling through thejet, it will perturb momentarily the stationnary shock structure. T ) In the toroidal case ( T ) (Fig. 5, top right), the interaction of themoving shock wave with the standing shocks has some simi-larities with the previous case. In the first rarefaction zone, theejecta accelerates adiabatically before it starts interacting withthe first shock. The resulting moving shock wave is subject to astrong radial expansion after this first interaction and slows down close to its initial value of γ ms =
10 at Z ∼ R jet . Then themoving shock sees its velocity increase between each recollima-tion shock (Fig. 6, green curve). Due to its initial higher inertiacompared to the surrounding jet, the ejecta undergoes a strongerinteraction when its crosses the recollimation zones. Therefore,the thermal pressure of the ejecta increases more than the ambi-ent flow. Afterwards, in the rarefaction zone, the higher pressureand inertia of the shocked ejecta starts to behave as a fireballwithin the surrounding jet. As a result, the moving shock waveinduced by the ejecta reaches a higher Lorentz factor than thesurrounding jet.As mentioned before (Section 4.2), the rarefaction zones arelarger than in the other cases, especially after the fourth standingshock, where the shock wave is strongly accelerated to a valueof γ ms (cid:39)
30 at Z ∼ R jet . Then the moving shock interactsstrongly with all the following standing shocks and causes themto oscillate along the jet axis (with a typical oscillation time closeto ∼ R jet ). P ) In the poloidal case ( P ) (Fig. 5, bottom left), the moving shockwave undergoes a strong acceleration to reach γ ms (cid:39)
21 beforeit interacts with the first stationary features and slows down to γ ms (cid:39)
13 at Z ∼ R jet (Fig. 6, orange curve.) In a secondphase, the resulting moving shock propagates through the suc-cessive rarefaction zones where it accelerates and the compres-sion zones where it decelerates. The mean Lorentz factor of themoving shock increases with distance until it reaches a station-ary shock at the distance Z ∼ R jet . This acceleration resultsfrom the expansion of the inner jet component in this region. Be-yond, the moving shock continues the propagation in the innerjet component which is compressed by the outer jet, and trans-verse instabilities grow along the stationary shock wave. In thisregion, the rarefaction zones are smaller and they are subject tomultiple small scale stationary shocks. As a result, the movingshock wave decelerates to reach a Lorentz factor of 14. As in thetoroidal case, beyond the fourth stationary shock, the passage ofthe moving shock induces an oscillation. HL ) In the helical case ( HL ) (Fig. 5, bottom right), the moving shock,after crossing the first internal shock wave at z = R jet , alsoundergoes a strong acceleration in the first rarefaction. Like inthe poloidal case, this strong acceleration is followed by a stronginteraction of the moving shock with stationary shocks, leadingto the deceleration of the moving shock to ( γ ms =
10) (Fig. 6,red curve.) Beyond the fourth standing shock, the moving shockaccelerates again to γ ms (cid:39)
6. Modeling the radiative transfer
In a post-processing step using a separate code, we evaluate thesynchrotron radiation of an electron population in the radio bandand solve the radiative transfer equation along a given line ofsight. We construct synchrotron emission maps to study the vari-ation of the flux observed in di ff erent zones along the jet overtime, as well as light curves of the spectral flux density integratedover the full simulated jet. In each cell, the relativistic electrons Article number, page 8 of 19. Fichet de Clairfontaine, Z. Meliani, A. Zech, O.Hervet: Flux variability from ejecta in structured relativistic jets. F i g . : S n a p s ho t: d i ff e r e n tt yp e s o f s t r u c t u r e d j e t s ( H , T , P a nd H L ) w it h a n i n j ec t e dp e r t u r b a ti on . T h e g e n e r a t e d m ov i ng s ho c k w a v e i s l o ca t e d a t ∼ R j e t fr o m t h e b a s e . T h e p r op a g a ti ono f t h e j e ti s go i ng fr o m l e f tt o r i gh t a l ong t h e i n c r ea s i ngv a l u e o f Z . F o r eac h ca s e , t h e d e n s it y c on t ou r( i n l og - s ca l e ) i s d r a w non t h e bo tt o m s i d ea nd t h e t h e r m a l e n e r gy c on t ou r on t h e t op s i d e . U n it s onx a ndy - a x i s i n R j e t un it . Article number, page 9 of 19 & A proofs: manuscript no. 2021_01_flare
Fig. 6: Evolution of the Lorentz factor of the moving shock as afunction of time in the co-moving frame. We represent, in di ff er-ent colors, the hydrodynamic ( H ), toroidal ( T ), poloidal ( P ) andthe helical case ( HL ) of jet.population is set with a power law, as expected for shock accel-eration, N e ( γ ) d γ = K γ − p d γ , (15)where γ min < γ < γ max . Following (Gómez et al. 1995), wedefine the normalization coe ffi cient as, K = e th , e (cid:0) p − (cid:1) − C − pE p − − C − pE n e (cid:0) p − (cid:1) p − , (16)where e th , e = (cid:15) e e th is the fraction of thermal energy carried bythe electrons, n e = (cid:15) e n is the fraction of the electron numberdensity with (cid:15) e = .
1, p = . ffi cient C E = γ max /γ min is set to 10 . As a simplification,we do not take into account radiative losses of the relativisticelectrons and assume, γ min = e th , e n e p − − − C − pE − C − pE . (17)In the present application, we are focusing only on the radioband, where the e ff ect of radiative cooling is the smallest.The specific intensity for each cell with index " i " is deter-mined in the frame of a stationary observer at the location of thesource (quantities in the co-moving frame are noted x (cid:48) ), I ν ; i = I ν ; i − exp ( − τ ν ) + S ν ; i (cid:0) − exp ( − τ ν ) (cid:1) , (18)where τ ν is the optical depth due to synchrotron self-absorptionand S ν is the synchrotron source function.Fig. 7 shows schematically how the contribution of each cellalong the line of sight to the total specific intensity I ν is esti-mated. It should be noted that here we do not account for lightcrossing time e ff ects, which can lead to a superposition of thesignal from di ff erent emission regions along the jet due to the Fig. 7: Schematic representation of the resolution of the radiativetransfer equation. We sum the di ff erent contributions along theline of sight.relativistic movement of the jet and the final speed of signalpropagation. In the radio band, this e ff ect is expected to be ofminor importance (Chiaberge & Ghisellini 1999).The specific intensity I ν (cid:16) in (cid:104) erg s − cm − Hz − sterad − (cid:105)(cid:17) depends on the specific emission coe ffi cient j ν , the absorptioncoe ffi cient α ν and optical depth τ ν . They are transformed to theobserver frame at the location of the source for each cell, j ν = δ j ν (cid:48) , (19) α ν = δ − α ν (cid:48) , (20) τ ν = τ ν (cid:48) . (21)These transformations (Rybicki & Lightman 1979) dependon the Doppler factor δ = ( γ (1 − β cos ( θ obs ))) − with γ = / √ − v the bulk Lorentz factor of the material in the cell and θ obs is the angle between the direction of the jet axis and the lineof sight. The di ff erent quantities ( j ν , α ν and τ ν ) are estimatedin each cell following the approximations given by (Katarzy´nskiet al. 2001) which are appropriate for the cases studied here. Thesynchrotron flux in the observer frame on Earth is determinedby, F ν = S e d l (1 + z ) I ν , (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 corre-sponding to M 87 ( z = . . R jet . To be able to distinguishbetween the emission from the jet and from the ejecta, we adjust Article number, page 10 of 19. Fichet de Clairfontaine, Z. Meliani, A. Zech, O.Hervet: Flux variability from ejecta in structured relativistic jets. an asymmetric 2D Gaussian distribution to the ejecta at eachtime 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 thepost-processing step. In this case, we assume that the magneticenergy density is a fraction (cid:15) B of the thermal energy density, B turb = √ (cid:15) B (cid:15) the , (23)with (cid:15) B = . Synthetic images are build by integrating the radiative transferequation (Eq. 18) along a line of sight, for all four cases ( H , T , P , HL ) and for one viewing angle θ obs = ◦ (Fig. 8). We show thelight curves realized for the four cases tested at θ obs = ◦ (Fig.9) and for the helical case at θ obs = ◦ , ◦ (Fig. 10). The lightcurve is computed by integrating the flux on the full simulationbox for the three observation angles. For an observation angle of θ obs = ◦ , it is also computed by integrating only the emissionissued from the moving shock wave using the method describedin Section 6.1.The synthetic image of the stationary hydrodynamic jet ( H )emission obtained for the observation angle θ obs = ◦ is pre-sented in Fig. 8 (top). The emission is dominated by the station-ary knots within the inner jet component. The intensity and thesize of the emitting knots decrease with distance as a result ofthe slow damping of shock strength with distance (Section 3).The contribution from the external jet component to the emis-sion is negligible, since the stationary shock wave within this jetcomponent is weak.The interaction between the moving shock wave, inducedby the ejecta, with the stationary shocks rapidly increases theinjection of accelerated particles and the associated synchrotronemission, causing a flare. The emission from the moving shockwave increases each time it passes through a standing knot.Moreover, after each collision, the stationary compression wavewithin the jet cannot ensure the radial cohesion of the shockedknot, thus the heated knot undergoes adiabatic expansion untilit reaches the interface between the inner-outer jet. Afterwardsthe knot cools down and contracts. A remnant emission isassociated to this adiabatic expansion. This illumination isvisible on the jet light curve obtained with observation angle θ obs = ◦ (Fig. 9, top left). It contributes to the emissionduring the decay phase of the observed flares (Fig. 9, top left).Moreover, each time a knot is shocked by the moving shockwave, it starts to stretch and to oscillate along the jet axis. Thisbehavior is associated with local changes of jet characteristics,such as the Mach number M and the inner jet radius R jet , in . Asa result, the structure of the standing shock waves evolves overa long period. This variability is associated with an increase ofthe base emission of the jet compared to the stationary case.The synthetic image of the stationary state for the toroidalcase ( T ) is shown in the second panel of Fig. 8. Like in the hydro-dynamic case ( H ), the emission comes mainly from the standingshocks in the inner jet, but the emission from the region betweenthe shocks is stronger. Moreover, the stationary shocks appearmore elongated along the jet and especially beyond the fourthknot. This behavior results from the toroidal magnetic field (Sec-tion 4.2). The moving shock within the jet increases the emission, andthe interaction with the stationary knots produces flares (Fig. 9,bottom left). In comparison with the hydrodynamic case ( H ), thestrength of the flares is variable with time. The strongest flare oc-curs at a time t ∼ R jet /δ when the moving shock crosses thefourth (and strongest) standing shock. This interaction leads toa strong deformation and oscillation of the knots along the jetaxis. The following relaxation of the knots is associated with aslow flux decline in the light curve and produces an asymmetricflare. This asymmetry occurs in the strongest flares and is en-hanced when a knot splits in two after interacting strongly withthe moving shock.The synthetic image for the poloidal case ( P ) is shown in thepanel on Fig. 8, where the emission of the stationary shock struc-ture is represented. The emission comes mainly from the part ofthe stationary shocks very near to the jet axis, where the shock issu ffi ciently strong. In comparison to cases H and T , the steadyshock wave is weaker in the poloidal case. A complex emissionstructure is visible due to the magnetic pressure along the Z -axisand the expansion of the outer toward the inner jet (Section 4.3).In this case, the variations in flux observed between shock andrarefaction zones are weak. Due to the expansion of the outerjet, the emission structure is more elongated along the jet axisand therefore less pronounced. As before in the P case, the mov-ing shock induced by the ejecta triggers flares as it crosses thesteady knots (Fig. 9, top right) but the e ffi ciency of its interac-tion is weaker than in the H and T cases. Indeed, the poloidalmagnetic field tends to damp the strength of the steady shockwaves and it limits the transverse expansion of the moving shockwave. The resulting light curve is characterized by weak flares.Furthermore, when the moving shock interacts with a standingshock, instabilities form along the steady shock wave. These in-stabilities will propagate through the jet in the form of severalmoving shocks and thus locally increase slightly the synchrotronemission.Finally, the synthetic images of the helical case ( HL ) canbe see in the last panel of Fig. 8, where the emission of thestationary shock structure is represented. As in all the otherscases, the emission comes mainly from the shock zones withinthe inner jet. However, the observed synchrotron flux is twice ashigh as in the other cases, due to the combination of the toroidalfield e ff ects with strong emission emanating from the center ofthe knots and the poloidal field e ff ects with di ff use emissionresulting from the compression of the inner jet by the outercomponent.As in the previous cases, the moving shock will increase locallythe synchrotron flux emitted and in particular when it crossesa stationary knot. But despite the fact that a toroidal magneticfield component is present, the shape of the light curve is verysimilar to the poloidal case ( P ) at θ obs = ◦ (Fig. 9, bottomright). As the emission is extended along the axis in the internaljet, variability is observed essentially from the inner jet itselfdue to the propagation of the moving shock. The presence of apoloidal field allows the development of transverse instabilitiesalong steady shock waves. These instabilities will propagate inthe inner jet in the form of moving shocks behind the principalmoving shock and increase the emission from the jet. Wenotice that the impact of a magnetic tension due to the toroidalcomponent is visible when the moving shock interacts withthe last shock zones. In fact, the moving shock interacts morestrongly with the last standing knot where the most pronouncedflare occurs. As before, the interaction of the moving shock withstanding knot induces an knot oscillation along the jet axis. Article number, page 11 of 19 & A proofs: manuscript no. 2021_01_flare
Fig. 8: Snapshot : synchrotron emission map of the di ff erent types of jets ( H , T , P , and HL ) without ejecta and stationnary. Eachmap represent the flux intensity in mJy unit. The R and Z -axis are given in R jet unit. These maps represent a resolution of theradiative transfer equation with an angle between the jet propagation axis and the line of sight equal to θ obs = ◦ and ν = H ) (Fig. 10) for θ obs = ◦ and θ obs = ◦ for 1 GHz. Compared to the one obtained for θ obs = ◦ , theyshow the expected increase in the overall flux due to strongerDoppler beaming. The flares are shorter in duration due to theDoppler e ff ect, but also due to self absorption. In certain cases,the dense ejecta partially hides the shocked knot behind it. Foran observation angle θ obs = ◦ , the self absorption e ff ect is verysignificant. The mean intensity decreases with distance, since,as the ejecta propagates within the jet, it hides a large number ofstationary knots.
7. Comparison with radio observations, the case of3C 273
Our study shows the complex link between ejecta propagatingin magnetized two-flow jets with stationary shocks, and its ob-served radio variability. All simulations in the present study con-sider a kinetic power of the outer jet larger than the one of theinner 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 cri-terion discriminating between types of VLBI radio jets, whichare themselves associated with spectral types of blazars (Hervetet al. 2016; Piner & Edwards 2018; Lister et al. 2019). Two-flowjets with kinetic powers within the same order of magnitude,such as the ones simulated in this study, were found to be themost similar to FSRQ-type blazars (flat spectrum radio quasar).In order to compare the results of our simulations withan astrophysical case, we focus on the radio observations ofone of the brightest and best monitored FSRQs over decades,3C 273 (B1226 + z = . / mas considering H =
70 km.s − .Mpc − . 3C 273 displays a peculiar mix of fastmoving and quasi-stationary radio knots. This hybrid radio kine-matic behavior is most often observed in intermediate blazars(Low or Intermediate frequency-peaked BL Lacs) (Hervet et al.2016). However 3C 273 significantly di ff ers from these sourcesin the way that quasi-stationary knots are visible only duringlow-activity periods of the jet, not continuously.For this study, we used 15.3 GHz observations from theVLBA, analyzed by the MOJAVE team up to August 2019 (Lis- Article number, page 12 of 19. Fichet de Clairfontaine, Z. Meliani, A. Zech, O.Hervet: Flux variability from ejecta in structured relativistic jets.
Fig. 9: Light curve obtained by integrating the total synchrotron flux emitted from a simulation box with a size of[R = , Z = R jet . The computation of the light curve is realized from the four di ff erent cases of jets ( H , T , P and HL ). The fluxis integrated from t = t ∼ R jet /δ with δ ( θ obs = ◦ ) = . ν = Article number, page 13 of 19 & A proofs: manuscript no. 2021_01_flare
Fig. 10: Light curve obtain by integrating the total syn-chrotron flux emitted from a simulation box with a size of[R = , Z = R jet . The computation of the light curve is real-ized for the hydrodynamic case ( H ) of jet. The flux is integratedfrom t = t ∼ (90 , R jet /δ for δ ( θ obs = ◦ ) (cid:39) .
57 (top) and δ ( θ obs = ◦ ) (cid:39)
18 (bottom) in the absolute frame and ν = . ± . a in Daly & Marscher (1988)). However study-ing the variability at higher energy, outside the synchrotron-self-absorbed frequencies, would be necessary to confirm this as-sumption.
8. Discussion
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 findthe well known diamond structure of standing shocks with aclear succession of shock and rarefaction zones in an over-pressured jet. A regular standing-shock structure with a linearintensity decrease over the distance from the base is observed.We found a combined e ff ect of the large-scale magnetic config-uration and jet stratification on the details of the shock struc-ture. As observed in Hervet et al. (2017), the jet stratification in-duces a larger variety of standing-shock structures due to the in-terferences between the characteristics waves from the two lay-ers of the jet. The toroidal magnetic field strengthens the stand-ing shocks by inducing intense rarefaction regions, where theplasma is strongly accelerated (Fig. 1 (top right)). Stratificationleads to the development of Rayleigh-Taylor instabilities at theouter-inner jet interface and along the standing shock wave asobserved in Toma et al. (2017) in the hydrodynamic case. In thepoloidal and helical cases, the magnetic pressure along the Z -axis amplifies these instabilities along the jet. They grow andheat both jet components. Within the inner jet, these instabilitiesinterfere with the standing shocks and lead to a smoother struc-ture and the appearance of a turbulent region at large distance(see Fig. 5 (bottom)).In a transverse structured jet, the presence of a structuredmagnetic field amplifies the jet opening angle compared to thehydrodynamic case (Fig. 3). This is especially apparent in thepoloidal and in the helical cases. Even if an outer jet layer shieldsthe inner part from a homogeneous ambient medium (Porth &Komissarov 2015), the magnetic field modifies the topology ofthe characteristic waves of the fluid. The presence of a toroidalmagnetic field component tends to limit this transverse expan-sion as observed in Mizuno et al. (2015) and at large distancesin our simulations. On the other hand, the poloidal magneticfield component induces instabilities as we saw; these instabili-ties lead to a jet decollimation at medium and large distances.It was shown by Martí (2015) that introducing an azimuthalvelocity component of the jet flow leads to a centrifugal forcethat further increases the jet opening. This e ff ect is not currentlytreated in our models, but we expect that it may have an impacton the standing shock structure, which we will evaluate in a fu-ture study. In all cases, the moving shock interacts with the successivestanding rarefaction zones, where it is accelerated adiabatically,
Article number, page 14 of 19. Fichet de Clairfontaine, Z. Meliani, A. Zech, O.Hervet: Flux variability from ejecta in structured relativistic jets. F l u x [ J y ] OVRO
MJD F l u x [ J y ] corek22k31k37k36k35 m a s t o t h e c o r e k22 k31 k37k39 k36k32 k35 Year
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 firmlyidentified 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 grayband 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 theobserved core variability is actually due to the flux increase of the emerging moving knots when they are indiscernible from the coredue to the limits of the VLBA angular resolution. Vertical lines show the most likely time when the moving knots pass through thestationary zone defined by k35 and the purple dashed line, with its associated uncertainty in gray.and collides with standing shocks, where it is heated and decel-erated. These interactions lead to the appearance of flares in thelight curves, due to thermal energy increase (cf. also Gómezet al. 1997; Fromm et al. 2016). The presence of a toroidal com-ponent of the magnetic field ensures the cohesion of the mov-ing shock and of the standing shocks. The fact that the movingand standing shock zones are very compact is reflected in theemission of intense and clearly marked flares. We recover a sim-ilar behavior in the hydrodynamic case. On the contrary, where a poloidal field is present, the interaction between the movingshock wave and the more di ff use steady shock zones is weakerand its associated flares are less pronounced.As we saw, the large variety of flares obtained is related to theintensity of the di ff erent knots, stemming from the combinationof the characteristic waves of the plasma. The outer jet compo-nent allows interferences between the stationary shock waves inthe inner jet and those of the outer jet. For certain conditions, thetwo stationary shock waves can combine and lead to a particu- Article number, page 15 of 19 & A proofs: manuscript no. 2021_01_flare
Fig. 12: Light curve obtain by integrating the total synchrotronflux emitted during the first three interactions between the mov-ing shock and the standing shocks. The computation of the lightcurve is realized for the four di ff erent cases of jets ( H , T , P and HL ) for δ ( θ obs = ◦ ) (cid:39)
18 and ν = . P and HL ) for δ ( θ obs = ◦ ) (cid:39)
10 and ν = . ff ect leads to astrong standing shock, with an important rarefaction zone behindit (cf. the fourth standing shock in our simulations). This stand-ing shock region is linked to the luminous flare emission in thelight curve (Fig. 9). After this interaction, the ejecta is stronglyaccelerated in the rarefaction zone and will lose its cohesion dueto fast adiabatic expansion. The simulated flares in H , T and P cases are characterized by atemporal 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 resolu-tion. 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 recollimationshocks, already observed by Agudo et al. (2001) and Mimicaet al. (2009), which will perturb standing knots along the jet axisand make them oscillate. This is associated with remnant emis-sion from the shocked knots during their adiabatic cooling phase.This process is observed in all cases, but most clearly in the caseswith pronounced interaction between the moving shock and thestanding shock (cases H and T ). Indeed, after the strong flare, wesee in the toroidal case a delay between the emission of the ejectaand the emission of the jet (Fig. 9). This is the emission signa-ture of the shocked knot, which causes a slight asymmetry. Thisadditional radiation counterpart will tend to soften the slope af-ter the ejecta has passed. This behavior is in accordance with theobserved flare structure described by Hovatta et al. (2008); Niep-pola et al. (2009) and obtained numerically from over-pressuredjets by Gómez et al. (1997); Fromm et al. (2016). As the angle decreases to 15 ◦ or 2 ◦ (Fig. 10), the e ff ect of theabsorption of the moving shock along the line of sight becomesmore important at low frequencies. Thus, the observed flare in-tensity, duration and asymmetry decrease as the occultation bythe moving shock becomes more important. Moreover, after theflare associated with the interaction of the moving shock withthe most intense steady knot, the flux intensity decreases. This isalso observed by Fromm et al. (2016) but, in our case, we havenot taken into account the e ff ect of the "light travel delay" ofphotons emitted in di ff erent parts of the jet, which would lead tosmoother light curves with a longer decay. This behavior is welldistinguishable with viewing angle θ obs = ◦ . With our shock propagation model, we arrive at a consistent, fornow qualitative interpretation of the flaring behavior observed in3C 273. Detailed modeling of the observed light curve is out ofthe scope of this work and will be addressed in a future publi-cation. The number of recollimation shocks present in our sim-ulated jets is much higher than the number of standing knotsobserved in the case of 3C 273, mostly due to the purely cylin-drical shape of our jets. In real jets that are imperfectly aligned,the superposition of knots along the line of sight and radio opac-ity might also lead to an obscuration of standing knots for theobserver. 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 maybe interpreted as series of recollimation shocks. Before the mov-ing shock waves interact with the first standing shock, extrapo-lation of the emitted flux seems to indicate that they are under-going a strong acceleration. This phenomenon appears clearlyin our results, with the moving shock passing through the firstrarefaction zone before interacting with the first knot. The accel-eration seems especially pronounced in the magnetized cases asshown in the Fig. 6. After the passage of moving shock waves,the disruption of the standing knots can be explained by the dy-namics and / or by the radiative transfer within the jet. Concern-ing the dynamics aspect, the trailing recollimation shocks canperturb momentarily the standing-shock structure, as we saw inour simulations. On the other hand, the apparent "disappearance" Article number, page 16 of 19. Fichet de Clairfontaine, Z. Meliani, A. Zech, O.Hervet: Flux variability from ejecta in structured relativistic jets. of standing knots, if confirmed, may suggest that they are simplyobscured by the brighter moving knots (k31 or k37) that hide thequasi-stationary knot, like we see in our synthetic light curve.In our model, the light curve (Fig. 12) obtained at frequency15.3 GHz (OVRO / MOJAVE frequency) and viewing angle 2 ◦ (e.g., Hervet et al. 2016) shows the interaction of the movingshock with the first two stationary shocks. The light curveintegrates the flux emitted by the whole jet like a single dishtelescope. We can thus compare our results with the OVROdata, in particular with the k37 event, which is isolated fromother events.We should note that the observations from Jorstad et al. (2017)found a di ff erent viewing angle of 6 . ± . ◦ . To investigate thee ff ects of the viewing angle, which is not precisely known, wealso compute the light curves for 6 ◦ (Fig.13).In the simulations, as in the observations, we find a flareduring each interaction between the moving shock wave anda standing shock. However, the typical flare duration is verydi ff erent ( ∼
800 days for k37 and ∼
150 days for our flares onaverage). This could be due to di ff erences in the morphologyof the jet, the size of the ejecta and stationary knots, and theuncertain value of the observation angle of 3C 273.The shape of the observed flares seems to show a small asym-metry. To quantify this e ff ect, we used the method proposed byRoy et al. (2018); Nalewajko (2013) to compare the doubling(or halving) time in the rise (or fall) phase of the flare. Applyingthe method to the k37 flare, we found ξ k37 = . ± . ξ H = . ± . ξ T = . ± . ξ P = . ± .
03 and ξ HL = − . ± .
06 for the hydrodynamic, toroidal, poloidal andhelical case for θ obs = ◦ (Fig. 12).At θ obs = ◦ , the shape of the light curve changes signifi-cantly due to beaming e ff ects, such that the partial superpositionof the first three flares does not allow us to clearly determinethe presence or absence of asymmetries. The peak of the secondflare is only well visible in the P and HL cases (Fig. 13). Inthe HL case, no asymmetry is found, while in the P case it hasswitched signs compared to the simulation at θ obs = ◦ (wefound for this case ξ P = − . ± . ff ects, which may play an important role at small angles.A more detailed study is beyond the scope of this work and willbe applied in a future study to dedicated simulations for a givendata set.The amplitude of the flares is on average much larger in3C 273 where the variability in radio can reach around twicethe baseline value, compared to an increase of ∼
15% in thesame-frequency, small angle simulations (Fig. 12). This di ff er-ence depends on the same characteristics as the flare duration.For the general picture, the main di ff erence with 3C 273 is thepresence of only two visible recollimation shocks, where onlythe first one is linked to strong flares, compared to a greaternumber of standing shocks and resulting flares in our simula-tions. The number of knots observed is however strongly linkedto the angular resolution and sensitivity of the VLBA. Anotherstationary knot very close to the core at ∼ . − .
16 mas, notedA1, ST1, or S1 has been detected when observing the jet at 43.2GHz (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 e ff ectmay be explained in 3C 273 by a strong interaction betweenshock waves in the inner jet and outer jet occurring at the posi-tion of k35, at 0.36 mas to the core. This zone can lead to a majoroutburst when interacting with a moving shock. It also has thespecificity of disrupting the downstream shock structure, whichwould explain the weak presence of the downstream stationaryshock k32, as well as the absence of significant flare event asso-ciated with it, and the apparent adiabatic cooling of ejecta mov-ing outside of k35.To directly model observed flares with our radiative SRMHDcode, a more important opening of the jet, reflecting changes inthe density profile of the ambient medium, will be required. Thisshould lead to a shock structure dominated by a few standingshocks close to the core. An implementation of the "light traveldelay" e ff ect (Chiaberge & Ghisellini 1999, cf.) and of radiativecooling will also be needed for a more realistic description of theradiative emission.
9. Conclusions
We have investigated the e ff ect of the large-scale magneticfield on the standing shocks and their interaction with ejectawithin a two-component relativistic jet. The associated lightcurves, which were computed at two radio frequencies ( ν = (1 , .
3) GHz) and for several observation angles ( θ obs = (2 ◦ , ◦ , ◦ )), show a variety of flares with varying durationsand amplitudes.Two-component magnetized jets are characterized by a com-plex standing-shock structure due to the interaction of character-istic waves propagating in the two jet components. In this way,jet stratification leads to the appearance of radio knots with arange of intensities along the jet. This is especially apparent inthe toroidal case, where we recover a strong standing shock, giv-ing rise to a pronounced flare in the interaction with the mov-ing shock wave. Temporal asymmetry associated to the relax-ation phase of the shocked standing knot is well visible for thestrongest flares. The introduction of large-scale magnetic fieldsis seen to cause an intrinsic opening of the jet with an openingangle up to three times larger than the hydrodynamic case forour jet configuration.Our scenario of moving ejecta interacting with standingshocks inside a two-component jet provides a good descrip-tion of the kinematics and light curves seen in the jet of theFSRQ type blazar 3C 273 with VLBI and single-dish radio ob-servations. In a preliminary study, at observation angle of 2 ◦ ,an asymmetry in the simulated flare profiles, with a fall-timethat is longer than the rise-time, was seen for the hydrodynamic,toroidal and poloidal cases, consistent with what is observed inthe OVRO data for this source. Acknowledgements.
The authors thank the anonymous reviewer for his peer re-view. The computations of the SRMHD results were carried out on the OC-CIGEN cluster at CINES in Montpellier (project named lut6216 , allocation A0050406842 ) and from the MesoPSL Cluster at PSL University in the Obser-vatory of Paris. This research has made use of data from the MOJAVE databasethat is maintained by the MOJAVE team (Lister et al. 2018), and from theOVRO 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 au-thors wish to thank M. Lister for providing preliminary MOJAVE data points forthis study. We wish also thank M. Lemoine for providing insightful discussionduring the study. Article number, page 17 of 19 & A proofs: manuscript no. 2021_01_flare
References
Agudo, I., Gómez, J.-L., Martí, J.-M., et al. 2001, The Astrophysical Journal,549, L183Asada, K., Nakamura, M., Doi, A., Nagai, H., & Inoue, M. 2014, ApJ, 781, L2Avachat, S. S., Perlman, E. S., Adams, S. C., et al. 2016, The AstrophysicalJournal, 832, 3Biretta, J. A., Junor, W., & Livio, M. 2002, New A Rev., 46, 239Blandford, R., Yuan, Y., Hoshino, M., & Sironi, L. 2017, Space Sci. Rev., 207,291Blandford, R. D. & Koenigl, A. 1979, Astrophys. Lett., 20, 15Blandford, R. D. & Königl, A. 1979, The Astrophysical Journal, 232, 34Blandford, R. D. & Payne, D. G. 1982, MNRAS, 199, 883Blandford, R. D. & Znajek, R. L. 1977, MNRAS, 179, 433Britzen, S., Kudryavtseva, N. A., Witzel, A., et al. 2010, A&A, 511, A57Cada, M. & Torrilhon, M. 2009, J. Comput. Physics, 228, 4118Chiaberge, M. & Ghisellini, G. 1999, MNRAS, 306, 551Daly, R. A. & Marscher, A. P. 1988, ApJ, 334, 539Dedner, A., Kemm, F., Kröner, D., et al. 2002, Journal of Computational Physics,175, 645Fendt, C., Porth, O., & Vaidya, B. 2012, Journal of Physics: Conference Series,372, 012011Fromm, C. M., Perucho, M., Mimica, P., & Ros, E. 2016, A&A, 588, A101Fromm, C. M., Perucho, M., Porth, O., et al. 2018, A&A, 609, A80Fuentes, A., Gómez, J. L., Martí, J. M., & Perucho, M. 2018, The AstrophysicalJournal, 860, 121Gabuzda, D. C., Reichstein, A. R., & O’Neill, E. L. 2014, MNRAS, 444, 172Ghisellini, G., Tavecchio, F., & Chiaberge, M. 2005, A&A, 432, 401Ghisellini, G., Tavecchio, F., Maraschi, L., Celotti, A., & Sbarrato, T. 2014, Na-ture, 515, 376Giroletti, M., Giovannini, G., Feretti, L., et al. 2004, The Astrophysical Journal,600, 127Gómez, J. L., Martí, J. M., Marscher, A. P., Ibáñez, J. M., & Alberdi, A. 1997,The Astrophysical Journall, 482, L33Gómez, J. L., Martí, J. M., Marscher, A. P., Ibáñez, J. M., & Marcaide, J. M.1995, The Astrophysical Journal, 449Hervet, O., Boisson, C., & Sol, H. 2016, A&A, 592, A22Hervet, O., Meliani, Z., Zech, A., et al. 2017, å, 606, A103Hovatta, T., Nieppola, E., Tornikoski, M., et al. 2008, å, 485, 51Jorstad, S. G., Marscher, A. P., Lister, M. L., et al. 2005, AJ, 130, 1418Jorstad, S. G., Marscher, A. P., Morozova, D. A., et al. 2017, ApJ, 846, 98Jorstad, S. G., Marscher, A. P., Smith, P. S., et al. 2013, ApJ, 773, 147Katarzy´nski, K., Sol, H., & Kus, A. 2001, å, 367, 809Keppens, R., Meliani, Z., van Marle, A., et al. 2012, Journal of ComputationalPhysics, 231, 718 , special Issue: Computational Plasma PhysicsKim, D.-W., Trippe, S., & Kravchenko, E. V. 2020, A&A, 636, A62Kim, J. Y., Krichbaum, T. P., Lu, R. S., et al. 2018, A&A, 616, A188Lemoine, M. & Pelletier, G. 2010, MNRAS, 402, 321Lemoine, M., Pelletier, G., Vanthieghem, A., & Gremillet, L. 2019a,Phys. Rev. E, 100, 033210Lemoine, M., Vanthieghem, A., Pelletier, G., & Gremillet, L. 2019b,Phys. Rev. E, 100, 033209Lisakov, M. M., Kovalev, Y. Y., Savolainen, T., Hovatta, T., & Kutkin, A. M.2017, MNRAS, 468, 4478Lister, M. L., Aller, M. F., Aller, H. D., et al. 2018, ApJS, 234, 12Lister, M. L., Homan, D. C., Hovatta, T., et al. 2019, ApJ, 874, 43Marscher, A. P. & Gear, W. K. 1985, The Astrophysical Journal, 298, 114Marscher, A. P., Jorstad, S. G., D’Arcangelo, F. D., et al. 2008, Nature, 452, 966Marshall, H. L., Miller, B. P., Davis, D. S., et al. 2002, ApJ, 564, 683Martí, J. M. & Müller, E. 2015, Living Reviews in Computational Astrophysics,1, 3Martí, J.-M. 2015, MNRAS, 452, 3106Martí, J. M., Perucho, M., Gómez, J. L., & Fuentes, A. 2018, International Jour-nal of Modern Physics D, 27, 1844011Mathews, W. G. 1971, ApJ, 165, 147Meliani, Z. & Keppens, R. 2007, å, 475, 785Meliani, Z. & Keppens, R. 2009, ApJ, 705, 1594Meliani, Z., Keppens, R., & Giacomazzo, B. 2008, A&A, 491, 321Meliani, Z., Sauty, C., Tsinganos, K., & Vlahakis, N. 2004, A&A, 425, 773Mignone, A. & Bodo, G. 2006, MNRAS, 368, 1040Mimica, P., Aloy, M. A., Agudo, I., et al. 2009, The Astrophysical Journal, 696,1142Mizuno, Y., Gómez, J. L., Nishikawa, K.-I., et al. 2015, The Astrophysical Jour-nal, 809, 38Nagai, H., Haga, T., Giovannini, G., et al. 2014, ApJ, 785, 53Nalewajko, K. 2013, Monthly Notices of the Royal Astronomical Society, 430,1324Nieppola, E., Hovatta, T., Tornikoski, M., et al. 2009, The Astronomical Journal,137, 5022 Pelletier, G., Gremillet, L., Vanthieghem, A., & Lemoine, M. 2019, Phys. Rev. E,100, 013205Perlman, E. S., Biretta, J. A., Zhou, F., Sparks, W. B., & Macchetto, F. D. 1999,ApJ, 117, 2185Piner, B. G. & Edwards, P. G. 2018, ApJ, 853, 68Plotnikov, I., Grassi, A., & Grech, M. 2018, MNRAS, 477, 5238Porth, O., Fendt, C., Meliani, Z., & Vaidya, B. 2011, The Astrophysical Journal,737, 42Porth, O. & Komissarov, S. S. 2015, MNRAS, 452, 1089Richards, J. L., Max-Moerbeck, W., Pavlidou, V., et al. 2011, ApJS, 194, 29Rieger, F. M. & Du ff y, P. 2019, ApJ, 886, L26Roy, N., Chatterjee, R., Joshi, M., & Ghosh, A. 2018, MNRAS, 482, 743Rybicki, G. B. & Lightman, A. P. 1979, Radiative processes in astrophysicsSavolainen, T., Wiik, K., Valtaoja, E., & Tornikoski, M. 2006, A&A, 446, 71Siemiginowska, A., Stawarz, Ł., Cheung, C. C., et al. 2007, The AstrophysicalJournal, 657, 145Sironi, L., Spitkovsky, A., & Arons, J. 2013, The Astrophysical Journal, 771, 54Sol, H., Pelletier, G., & Asseo, E. 1989, MNRAS, 237, 411Strauss, M. A., Huchra, J. P., Davis, M., et al. 1992, ApJS, 83, 29Tavecchio, F. 2020, arXiv e-prints, arXiv:2011.03264Tavecchio, F. & Ghisellini, G. 2008, Monthly Notices of the Royal AstronomicalSociety: Letters, 385, L98Tavecchio, F. & Ghisellini, G. 2014, MNRAS, 443, 1224Toma, K., Komissarov, S. S., & Porth, O. 2017, MNRAS, 472, 1253Valtaoja, E., Lahteenmaki, A., & Terasranta, H. 1992, Astronomy & Astro-physics, Suppl., 95, 73Walker, R. C., Hardee, P. E., Davies, F. B., Ly, C., & Junor, W. 2018, The Astro-physical Journal, 855, 128Wehrle, A. E., Grupe, D., Jorstad, S. G., et al. 2016, The Astrophysical Journal,816, 53Wehrle, A. E., Marscher, A. P., Jorstad, S. G., et al. 2012, ApJ, 758, 72Wilson, A. S. & Yang, Y. 2002, ApJ, 568, 133Zanotti, O., Rezzolla, L., Del Zanna, L., & Palenzuela, C. 2010, A&A, 523, A8 Article number, page 18 of 19. Fichet de Clairfontaine, Z. Meliani, A. Zech, O.Hervet: Flux variability from ejecta in structured relativistic jets.
Appendix A: Additional Figures
Fig. A.1:
Left : Initial representation of the toroidal magnetic field strength in space following Eq. 13 and assuming (cid:16) B ϕ, in , , B ϕ, out , (cid:17) ≡ Right : Initial representation of the magnetic pressure following Eq. 14.: Initial representation of the magnetic pressure following Eq. 14.