Polarization Signatures of Relativistic Magnetohydrodynamic Shocks in the Blazar Emission Region - I. Force-free Helical Magnetic Fields
aa r X i v : . [ a s t r o - ph . H E ] D ec Submitted to
The Astrophysical Journal
Polarization Signatures of Relativistic MagnetohydrodynamicShocks in the Blazar Emission Region - I. Force-free HelicalMagnetic Fields
Haocheng Zhang , , Wei Deng , , Hui Li and Markus B¨ottcher ABSTRACT
The optical radiation and polarization signatures in blazars are known to behighly variable during flaring activities. It is frequently argued that shocks arethe main driver of the flaring events. However, the spectral variability modelingsgenerally lack detailed considerations of the self-consistent magnetic field evolu-tion modeling, thus so far the associated optical polarization signatures are poorlyunderstood. We present the first simultaneous modeling of the optical radiationand polarization signatures based on 3D magnetohydrodynamic simulations ofrelativistic shocks in the blazar emission environment, with the simplest physicalassumptions. By comparing the results with observations, we find that shocksin a weakly magnetized environment will largely lead to significant changes inthe optical polarization signatures, which are seldom seen in observations. Hencean emission region with relatively strong magnetization is preferred. In suchan environment, slow shocks may produce minor flares with either erratic po-larization fluctuations or considerable polarization variations, depending on theparameters; fast shocks can produce major flares with smooth PA rotations. Inaddition, the magnetic fields in both cases are observed to actively revert to theoriginal topology after the shocks. All these features are consistent with obser-vations. Future observations of the radiation and polarization signatures willfurther constrain the flaring mechanism and the blazar emission environment.
Subject headings: galaxies: active — galaxies: jets — gamma-rays: galaxies —radiation mechanisms: non-thermal — relativistic processes — polarization Astrophysical Institute, Department of Physics and Astronomy,Ohio University, Athens, OH 45701, USA Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA Department of Physics and Astronomy, University of Nevada Las Vegas, Las Vegas, NV 89154, USA Centre for Space Research, North-West University, Potchefstroom, 2520, South Africa . Introduction
Blazars are the most violent active galactic nuclei. They are known to emit nonthermal-dominated radiation from radio to γ -rays with strong variability across the entire electro-magnetic spectrum. It is generally agreed that the emission comes from an unresolved regionin a relativistic jet that is directed close to our line of sight (LOS). The blazar spectrumhas two components. The low-energy component, from radio to optical/UV, is known to bepolarized, with the polarization percentage ranging from a few to tens of percent. This is inagreement with the synchrotron emission from nonthermal electrons in a partially orderedmagnetic field. Several papers have demonstrated that the observed polarization signa-tures may indicate a helical magnetic structure (Lyutikov et al. 2005; Pushkarev et al. 2005;Zhang et al. 2015). The high-energy component, from X-rays to γ -rays, is usually inter-preted as originating from inverse Compton scattering by the same nonthermal electrons ofsoft seed photons (e.g., Marscher & Gear 1985; Dermer et al. 1992; Sikora et al. 1994), but ahadronic origin cannot be ruled out (e.g., Mannheim & Biermann 1992; M¨ucke & Protheroe2001; B¨ottcher et al. 2013). Both spectral components of blazars exhibit fast variability.Many observations show flares in various observational bands lasting days or even hours(e.g., Ciprini 2011; Chatterjee et al. 2012); in particular, high-energy γ -rays sometimes showvariability time scale within several tens of minutes (e.g., Aharonian et al. 2007; Albert et al.2007). Based on the causality relation, R . δct , where δ is the Doppler factor and t is theobserved flaring time scale, the size of the emission region R in the comoving frame of theemission region should generally be within 0 . pc . Therefore, it is often argued in the blazarspectral variability fittings that the blazar emission region is a small region near the broadline region of the blazar jet (e.g., Tavecchio et al. 2010; B¨ottcher et al. 2013; Barnacka et al.2014).Blazar spectral fittings have shown that the nonthermal electron spectra responsi-ble for the common blazar SEDs require power-law indices of & ∼ ∼ & ◦ ) PA rotations (e.g., Marscher et al. 2008; Abdo et al. 2010;Chandra et al. 2015). These flare + PA rotation events often feature smooth PA rotationsand apparently time-symmetric light curves and PD patterns. In addition to the individualobservations, Blinov et al. (2015) are performing a polarization monitoring program of alarge sample, in which they find that blazars with detected PA rotations generally exhibitstronger variations in the PA, and the rotations can go in both directions. Moreover, theyclaim that the PA rotations are probably physically connected to the γ -ray flares.In terms of modeling, several mechanisms have been put forward to understand thepolarization signatures, especially the PA rotations. An emission region following a helicaltrajectory (e.g., Marscher et al. 2008) can give rise to a smooth PA rotation, but such amodel prefers all rotations present in the same blazar to follow the same rotating direction,which contradicts observations. An initially chaotic magnetic field structure compressedby a flat shock (Laing 1980) or multiple small shocks (Marscher 2014) can produce erraticpolarization patterns and occasionally large polarization variations through random walks,but the resulting patterns are normally whimsical, and it has been shown that the observedPA rotations are unlikely to result completely from random walks (Kiehlmann et al. 2013;Blinov et al. 2015). Chen et al. (2014) and Zhang et al. (2014) have proposed a model inwhich a disturbance propagates through a cylindrical emission region pervaded by a he-lical magnetic field. This model can give rise to systematic PA rotations and apparentlytime-symmetric light curves and PD patterns. Based on this model, Zhang et al. (2015)and Chandra et al. (2015) have presented simultaneous fittings of SEDs, light curves andpolarization signatures for two flare + PA rotation events. Nonetheless, in all the modelsmentioned above, the magnetic field evolution is treated in an ad hoc way. Furthermore,3one of them can so far explain the statistical properties of the polarization variations asshown in Blinov et al. (2015).In this paper, we present the first polarization-dependent radiation modeling basedon relativistic magnetohydrodynamics (RMHD) simulations of shocks in helical magneticfields in the blazar emission region. This new approach performs detailed analysis on theinteraction between the shock propagation and the magnetic field evolution, which enablesus to constrain the blazar emission environment and the shock parameters. We demonstratethat the blazar emission region largely possesses significant magnetic energy. In such anenvironment, fast shocks can be the driver of flaring events, producing both the typicalpolarization fluctuations and the occasional large PA rotations. Our simulation results areconsistent with the statistical polarization properties reported in Blinov et al. (2015). Wewill describe our model setup in Section 2, illustrate the interaction between the shock andthe magnetic field as well as the consequent radiation and polarization signatures in Section3, present additional parameter studies in Section 4, and discuss the results in Section 5.
2. Model Setup
The purpose of this paper is to study the time-dependent radiation and polarizationsignatures resulting from an RMHD shock in the blazar emission environment, with thesimplest physical assumptions. In this section, we will describe our physical assumptionsand corresponding code structure in detail.
The blazar emission region is often considered to be an unresolved region near the broadline region of the jet. Generally speaking, the magnetic field topology and evolution as well asplasma flow dynamics in the blazar emission region should be linked with those in the large-scale jet. However, the oberved fast variability and high luminosity suggest that the emissionregion is an extraordinary, very small and localized region with fast evolution. Therefore,we argue that within the time scale of an individual flare, the blazar emission region canbe considered as uncorrelated with the properties of the large-scale jet, even though it isspatially embedded within the jet. On the other hand, the emission region is expected toremain well localized in the jet during flares due to pressure provided by the large-scale jetstructure. We assume that flares are due to a disturbance propagating through the emissionregion. Many observations show that light curves and even the time-dependent polarization4ignatures appear symmetric in time (e.g., Abdo et al. 2010). This means that the emissionregion and the disturbance are likely to possess some kind of symmetry. Zhang et al. (2014,2015) have shown that a helically symmetric emission region and disturbance, with detailedconsideration of light-travel-time effects (LTTEs), can naturally explain the apparently time-symmetric light curves and time-dependent polarization signatures. Therefore, we continueto use that assumption here.Due to the relativistic aberration, even though we are observing blazars nearly along thejet in the observer’s frame (typically, θ obs, ∼ / Γ ), where θ obs, and Γ are the angle betweenthe LOS and the jet direction and the Lorentz factor of the emission region in the observer’sframe, respectively, the angle θ obs between LOS and the jet axis in the comoving frame islikely around 90 ◦ (if θ obs, ∼ / Γ , then θ obs ∼ ◦ ). Zhang et al. (2014) have shown that ifthe LOS is fixed at some other angles, the general trend of polarization variations is onlyweakly affected. In addition, as is suggested in the bending jet scenario of the polarizationsignatures, a change in the LOS direction may significantly affect the polarization signatures.However, those effects are beyond the scope of our first-step study. Thus in all of thefollowing simulations, we choose θ obs = 90 ◦ in the comoving frame, and hence the Dopplerfactor δ ≡ (Γ [1 − β Γ cos θ obs, ]) − ∼ Γ .We use ideal RMHD simulations to describe the evolution of the system. The underlyingmodel assumes that a plasma jet is traveling with a Lorentz factor of Γ in the observer’sframe. In the comoving frame of the jet, a flow of plasma which is pervaded by a helicalmagnetic field travels at a Lorentz factor of Γ. Inside this flow lies a cylindrical emissionregion containing nonthermal particles. Along the path of the flow it encounters a flatstationary layer of the plasma (the disturbance), forming a shock wave. Thus the totalLorentz factor of the emission region in the observer’s frame is Γ = Γ × Γ . Although in thefollowing sections our RMHD simulations use different Lorentz factors Γ of the disturbancein the comoving frame of the emission region, Γ of the emission region in the observer’sframe is not constrained. Therefore, we assume a fixed Lorentz factor of Γ = 20 for theemission region in the observer’s frame, which is commonly inferred from SED modelings.Before the emission region interacts with the disturbance, the helical magnetic field isassumed to be in force-balance. The advantage of this magnetic field setup is that it willnaturally give rise to comparable poloidal and toroidal contributions, resulting in a relativelylow background PD without the need of any turbulence. In addition, this can simplify oursetup for the plasma density and the thermal pressure, which are taken to be uniform insidethe emission region. In this way, we assume that initially all flow conditions and magneticfields are laminar, with no preexisting turbulent components. Also, as the RMHD requiresan input for the thermal pressure, we assume that in the beginning the plasma is cold and5ence the thermal pressure is very small compared to the kinetic energy or the magneticenergy.In the comoving frame of the emission region, the disturbance and the resulting shockwill propagate through the emission region, and change the physical conditions as well asaccelerate particles at its location, generating a flare. We apply the simple assumption thatthe shock will inject fresh nonthermal electrons at the shock front whose energy is equal toa fixed small amount of the local shock kinetic energy. In this way, any deceleration and de-formation of the shock due to the magnetic field obstruction can be taken into consideration.In Section 4 we will show that the injection rate does not affect the polarization signaturestoo much.We briefly summarize our model assumptions in the following:1. The emission region is a small cylindrical region embedded in the large-scale jet.2. The evolution of the emission region is detached from the large-scale jet.3. The boundary of the emission region is held by a pressure wall, probably provided by thelarge-scale jet.4. In the comoving frame of the emission region, the observer is observing from the side ofthe emission region ( θ obs = 90 o ).5. The initial magnetic field inside the emission region is a force-free helical magnetic field.6. Initially all flow conditions and magnetic fields are laminar.7. Initially the plasma is cold.8. Initially the disturbance is a flat cylindrical region traveling relativistically in the comovingframe of the emission region.9. The disturbance will generate a shock in the emission region. At the shock front, theshock will inject nonthermal particles whose total energy is a constant, small fraction of thelocal shock kinetic energy. The above model is realized by the combination of the 3D multi-zone RMHD codeLA-COMPASS developed by Li & Li (2003) and the 3D multi-zone polarization-dependentray-tracing code 3DPol developed by Zhang et al. (2014). Fig. 1 (left) shows a sketch of thesimulation setup, and Table 1 presents some major parameters. The RMHD simulation isperformed in the comoving frame of the emission region, using Cartesian coordinates. Weapply the outflow boundary conditions for the simulation domain; in the following, we willdemonstrate that our setup should not suffer from any boundary effects. The simulationdomain is pervaded by a helical magnetic field (Fig. 1 (left)). In reality, the poloidal6agnetic field lines have to be closed, but based on our assumption the returning magneticflux is generally outside the simulation domain. The helical magnetic field is assumed to bein force-balance, in the form of B z = B × J ( kr ) B φ = B × J ( kr ) (1)where J and J are the Bessel functions of the first kind, r is the radius, and B and k arenormalization factors of the magnetic field strength and the radius, respectively. The B r component is assumed to be zero. In order to avoid any boundary effects in x , y directionsand comply with the assumed generally monodirectional B z , we add an exponential cutofffor both components when B z approaches zero ( kr ∼ . Parameters Length Time VelocityRelation L L /c c Code Unit 1 1 1Physical Value 8 . × cm . × s × cm s − Parameters Magnetic Field Thermal Pressure Plasma DensityRelation B B / (4 π ) B / (4 πc )Code Unit 1 1 1Physical Value 0 . G . × − erg cm − . × − g cm − Bulk Lorentz factor Γ θ obs ( ◦ ) 90Electron minimal energy γ Electron maximal energy γ × Electron power-law index p Table 1: Summary of parameters. Top: Conversion between the RMHD code units and thephysical value. Bottom: Additional parameters used in the 3DPol simulation. All parametersin this table are in the comoving frame of the emission region, except the bulk Lorentz factor.The emission region is set to be a fixed volume in the simulation domain ( r ∼ . z ranges from − . .
5, see Fig. 2). Since the simulation domain is in the comoving frameof the emission region, the disturbance will travel at a Lorentz factor of Γ in the z direction.We assume that the disturbance is a thin layer (0 . − − .
5) toallow some time for the shock wave to form in the numerical simulation (see Fig. 1 (left) andFig. 2). In order to avoid boundary effects in z direction, we use a sufficiently large rangein z ( −
20 to 20), so that at the end of the simulation, no signal has reached the z boundary.We define a magnetization factor in the emission region, σ = E em h (2)7here E em = B + E π is the electromagnetic energy density, h = ρc + ˆ γp ˆ γ − is the specificenthalpy, ρ is the plasma density, ˆ γ is the adiabatic index and p is the thermal pressure. Wefurther assume that the plasma is cold. Hence we choose ˆ γp ˆ γ − to be a small part (1 /
30) ofthe electromagnetic energy initially. In this way the magnetization factor is approximately σ ∼ E em ρc . The relativistic Alfv´en speed can then be expressed as V A ∼ c p πρc /B + 1 (3)Since the entire simulation domain is axisymmetric, we present a cut in the xz plane inFig. 2 to illustrate the initial condition, including the magnetic field and the velocity of thedisturbance. Notice that the B x component is not shown, as in a helical setup, in the xz plane this component is trivial. Also this figure is for an initial Γ = 10; for other Γ thevelocity appears slightly different. As the RMHD simulation takes Cartesian coordinates,we will transform the results to the cylindrical coordinates and feed into the 3DPol code.The 3DPol code is focused on the polarization-dependent synchrotron radiation. Ituses a cylindrical geometry, and further divides it evenly into multiple zones in r , φ and z directions. Each zone has its own magnetic field evolution, which is provided by the RMHDsimulations. The nonthermal electrons in each zone are assumed to range between a minimaland a maximal Lorentz factor γ and γ , with a fixed spectrum of power-law index of − n ( γ ) = nγ − e − γγ , f or γ > γ (4)The nonthermal electron normalization factor n is assumed to have two components: a fixeduniform initial background density n b , and an injected density n i during the shock. As theelectron spectrum is fixed, we will only consider the optical light curves and polarizationpatterns, which do not suffer from any synchrotron-self absorption and Faraday rotationeffects for our parameters. As in the RMHD simulation, all calculations are performed inthe comoving frame of the emission region. Given the local magnetic field information andthe nonthermal electron population, the 3DPol code will calculate the Stokes parameters atvarious frequencies at every time step in each zone. Next the code will employ ray-tracing totake account of all LTTEs (see Fig. 1 (right) for an illustration), and add up the incoherentStokes parameters arriving at the observer at the same time, then Lorentz transform tothe observer’s frame so as to obtain the total time-dependent radiation and polarizationsignatures.We define the PA in the following way. When the electric vector is parallel to the emis-sion region propagation, P A = 0 (toroidal component dominating). Since the PA has 180 ◦ ambiguity, the toroidal dominance happens at P A = 2 N × ◦ , and the poloidal dominanceat P A = (2 N + 1) × ◦ , where N is an integer.8 . Interaction between Shock and Magnetic Field In this section, we will investigate how the shock speed and the magnetization in theblazar emission environment will affect the time-dependent radiation and polarization signa-tures. Since we generate the shock by a relativistic disturbance, and the simulation is donein the comoving frame of the emission region, the shock front and its motion can hardly bedistinguished from the disturbance. Moreover, as we will see in the following results, duringthe propagation, both the disturbance and the shock speeds may vary in time. Hence weuse the initial Lorentz factor Γ of the disturbance rather than the actual time-dependentshock speed for the parameter study. For the same reason, the magnetization factor σ ′ inthe shock frame will vary in time as well; we choose the initial σ in the emission regioninstead. This σ is also a direct indicator of how much the emission region is magnetized atthe beginning. By comparing the results with the general observational features, we will beable to constrain the physics in the emission environment. For this purpose, we will presentthree cases, namely, σ ∼ .
01, Γ = 10 (Case I), σ ∼
1, Γ = 10 (Case II) and σ ∼
1, Γ = 3(Case III). The RMHD simulation results are shown in Figs. 3, 5 and 6 for Case I, II andIII, respectively, along with the associated radiation and polarization signatures in Fig. 7.For all the cases studied, we fix the initial nonthermal electron density n b and the character-istic magnetic field strength B . The magnetization factor is adjusted through changing theplasma density. Since we employ a simple electron spectrum, we only study the light curvesand polarization patterns in the optical band as examples. Additionally, as the backgroundnonthermal electron density and the injection rate are set as input parameters, we use therelative flux level, where 1 is approximately 10 − erg cm − s − . We first examine a kinetic energy dominated environment, Case I. Such an environmentis frequently assumed in a shock-initiated flaring model. We first notice that in a triviallymagnetized setting ( σ ∼ . V A ∼ . c , is much slower thanthat of the shock/disturbance. This speed is even slower inside the disturbance. As a result,the entire simulation domain can be treated similar to a hydro setup. Consequently, we canexpect that the shock/disturbance will strongly compress the plasma at the shock front, andpush it along with the propagation. However, in the ideal MHD, the plasma and the field linesare coupled (frozen-in field lines). Therefore, the toroidal magnetic field lines at the shockfront will be greatly compressed, leading to a much larger toroidal contribution; meanwhile,far downstream of the shock, the toroidal lines will be sufficiently stretched, giving rise toa perfectly ordered poloidal magnetic field (see Fig. 3, B y ). The physical effects described9bove are the consequence of a shock in a nearly hydro environment with the assumptionof the frozen-in field lines. Therefore, we expect that the above shock modifications to themagnetic field structure is not limited to our initial magnetic field topology, but will happenin other magnetic field setups as well. Notice that, however, due to the shock compression,the initial force balance in the radial direction between the magnetic field and plasma pressureand density may be broken. This will lead to some hydrodynamical perturbations in theplasma density and pressure in the radial direction, which are explicitly illustrated in Fig.4. However, those perturbations are only transported at sonic speeds. This is supported byFig. 4 where any plasma velocities in the radial direction that arise from those perturbationsare non-relativistic. Therefore, although given enough time they may significantly modifythe post-shock structures, they are much slower than the time scale that we are interested inin our RMHD simulation. Therefore, we find that those hydrodynamical effects are of minorimportance to our current study, and will not discuss them in the following.Case II has a moderately magnetized environment ( σ ∼ r c = 6(1 + σ ′ )1 + 2 σ ′ + √ σ ′ + 16 σ ′ + 1 (5)where r c and σ ′ are the shock compression ratio and the magnetization factor in the shockframe, respectively. As a result, at the shock front, we can see from the simulation that theshock compression becomes weaker due to the higher magnetization (Fig. 5, B y ). Further-more, although initially the emission region has a negligible magnetic force, the enhancementin the toroidal component at the shock front will break down this balance, resulting in acontracting magnetic force in the r direction, given by F B = − ∂∂r B φ + B z − B φ r (6)This contraction is transported by Alfv´en waves, which are mildly relativistic in this case.Therefore, in addition to the compression in the z direction at the shock front, there exists acontraction of the plasma in the r direction. Because of the frozen-in field lines, the poloidalcomponent will be strengthened at the shock front (Fig. 5, B z ). We can see from thesimulation that at a later stage, the increase in the poloidal component gradually gets closerto that in the toroidal component, and even the shape of the shock/disturbance is modified,showing a bullet shape at the central part (Fig. 5 V z ). Moreover, at variance with Case I,far downstream of the shock/disturbance, even though some plasma is pushed away by theshock/disturbance, the magnetic field is observed acting to revert to its initial topology (Fig.5, B y and B z ). 10n the spectral variability fitting, it is frequently argued that during the shock propaga-tion, the local magnetic field either generally maintains its strength and topology (Sokolov et al.2004; Joshi & B¨ottcher 2007), or reverts to its initial state after the shock (Chen et al. 2014;Zhang et al. 2014). This requirement is often necessary to produce the best fittings. Oursimulations demonstrate that this can only happen in an adequately magnetized emission re-gion. Therefore, detailed analysis of the magnetic field evolution and shock dynamics shouldplay an essential part in those spectral variability models.We notice that in both simulations, there are regions with negative V z (Figs. 3 and 5, V z ). This is likely due to the reverse shock. In reality, the reverse shock may also lead tosome radiation and polarization signatures, but its effect is expected to be weaker than themain shock. Thus in the following we will not consider it.Now we will consider the radiation and polarization signatures resulting from the aboveRMHD simulations, calculated using the 3DPol code. The general trends are similar forboth cases. Before the shock moves in, the entire emission region is axisymmetric. Althoughthe poloidal and toroidal components are comparable in the force-free setup, due to theLOS effect, the projected toroidal component onto the plane of sky is weaker than theprojected poloidal component. Therefore initially the emission region has a total polarizationdominated by the poloidal component (Fig. 7). When the shock moves in, the toroidalcomponent is strongly increased and fresh nonthermal electrons are injected at the shockfront. However, due to the LTTEs, only a small elliptical region on the right that is nearthe observer is seen (the near side, see Fig. 1 (right), red). Therefore, the flux is seen togradually increase. Nevertheless, since the toroidal enhancement is strong, a small flaringregion is adequate to push the polarization to be dominated by the toroidal component. As aresult, the PA quickly rotates to 180 ◦ and forms a plateau, which corresponds to the toroidaldomination. The PD first experiences a drop which indicates a switch in the dominationbetween the two components, then climbs up to a high level, implying a strong toroidaldominance (Fig. 7). When the shock is about to leave the emission region, the flaring regionreaches maximum (Fig. 1 (right), green). Hence the flux peaks. After that, the flaringregion moves far from the observer to the left (the far side, Fig. 1 (right), blue), so theflux gradually decreases in an apparently symmetric pattern. However, as some plasma andthe coupled magnetic field lines have been pushed away by the shock/disturbance, the totalmagnetic field strength is smaller than the initial state. Hence the flux level is lower thanthe initial value. For the same reason, the toroidal contribution is weakened near the end ofthe flare, hence the PD rises up higher than the initial value, revealing a stronger poloidalcontribution. Since the helical magnetic field direction on the far side is opposite to thaton the near side, the PA instead completes a 180 ◦ rotation to the initial value (180 ◦ PAambiguity, Fig. 7). 11he major differences between Case I and II are the following. First, due to the strongermagnetic field, Case II has a weaker compression in the toroidal component at the shockfront than Case I, while it experiences a rise in the poloidal contribution. Therefore, thePA rotation in Case II appears smoother and has a shorter toroidal dominated plateau thanCase I. Also the maximal flare level in Case II is lower. More importantly, near the endof the flare, the moderate magnetization in Case II is acting to restore the initial magnetictopology, hence even if the PD rises above the initial value, it immediately starts to decrease,indicating the recovery of the toroidal component (Fig. 7). Abdo et al. (2010) have shownthat at the end of a flare + PA swing event in 3C 279, the PD rises up above the initialvalue, which is immediately followed by a significant restoration to approximately the initialPD level. The same restoring phases are seen in a number of flare + PA rotation events aswell (Larionov et al. 2008, 2013; Morozova et al. 2014). On the other hand, in Case I themagnetization is too weak, thus the PD rises up to ∼
70% at the end of the flare and showsno restoration, indicating a purely poloidal magnetic structure. As is mentioned above, sucha phenomenon is generally the consequence of a trivially magnetized environment, and it isexpected to happen with other initial magnetic field topologies. Therefore, we argue thatsuch drastic changes in the polarization signatures are unrealistic compared to observations,so that a nearly unmagnetized emission environment model can be ruled out.
The effects of the shock speed are rather straightforward. We will illustrate those byexamining Case III, which shares the same magnetization as Case II but has a slower dis-turbance. A slower disturbance has less kinetic energy. Therefore, as the shock/disturbancepropagates, it cannot exert a very strong pressure onto the plasma; in return the plasma willdecelerate the shock/disturbance. Thus we see in Fig. 6 that the shock/disturbance is con-siderably slower than the initial value (Γ = 3). Therefore, the flare duration is much longerthan the previous cases (Fig. 7). Also the shock/disturbance will provide less compressionin the plasma at the shock front. As a result, we observe in the simulation that the toroidalcomponent is only moderately increased and stretched at the shock front and downstream ofthe shock/disturbance, respectively (Fig. 6, B y ). Moreover, the shower shock/disturbancetakes longer time to propagate through the emission region, allowing for further informationexchange by Alfv´en waves. Therefore, the contraction of the plasma at the shock front, dueto the enhanced toroidal component, is seen to catch up with the shock compression. Thisleads to a considerably strengthened poloidal component at the shock front. Consequently,the PA rotation disappears; instead it shows fluctuations around the initial value. However,at the flare top, the poloidal contribution is only marginally stronger than the toroidal con-12ribution. Therefore, the PD still displays a decrease at the middle of the flare (Fig. 7).Meanwhile, downstream of the shock/disturbance, Alfv´en waves restore the initial magneticfield topology to a large extent (Fig. 6, B y and B z ). Hence the final PD is not too muchhigher than the initial value, and gradually recovers after the flare (Fig. 7). Another inter-esting effect is that due to the weaker kinetic energy in the shock/disturbance and longerinfluence of Alfv´en waves, the shape of the shock/disturbance is distorted (Fig. 6, V z ). Thusthe polarization patterns, especially the PA, appear less symmetric in time than the previouscases (Fig. 7).
4. Parameter Study
In this section, we will perform parameter studies to constrain the shock speed and themagnetization of the emission environment, by comparing the polarization signatures withthe general observational properties. Since the injection rate at the shock front is somewhatarbitrary, we will not compare the light curves. Although in principle the injection rate willaffect the polarization signatures, in the following we will demonstrate that the polarizationsignatures are mainly the consequence of the magnetic field evolution. Here we use moderate(Γ = 6) and slow (Γ = 3) speeds, as well as weak ( σ ∼ . σ ∼
10) and moderate( σ ∼
1) magnetization factors as examples. σ ∼ . σ ∼ . B z and V z ).As a consequence, the polarization patterns appear relatively smooth, in particular the PArotation (Fig. 9). However, downstream of the shock/disturbance, Alfv´en waves fail torestore the initial magnetic topology. Hence at the end of the flare, the PD rises up andmaintains a high level ( ∼ .2. Strongly Magnetized Environment, σ ∼ σ ∼ B y and B z ). Therefore, the PApatterns for both cases exhibit no rotation (Fig. 11). Additionally, the shock/disturbance isgreatly distorted. Hence both the PD and the PA patterns appear rather erratic, showinga lot of bumps throughout the flare (Fig. 11), even though the general geometry is rathersymmetric. Finally, the restoration downstream of the shock/disturbance is very fast (Fig.10 B y and B z ), thus the PD only rises up gently above the initial value, and quickly recoversat the end of the flare (Fig. 11).Notice that for a sufficiently fast shock/disturbance, a PA rotation may reappear. How-ever, we expect that such a rotation is likely to be very smooth, without a long toroidaldominated plateau as in Case I to III. Meanwhile the PD variation can be within somesmall value. This may agree with the observed flare + PA rotation events. To summarize, astrongly magnetized emission environment may generate the typical erratic perturbations inthe polarization signatures in the condition of a relatively slow disturbance, together withthe smooth flare + PA rotation events provided a relatively fast disturbance. Thus we preferan emission region with high magnetization. σ ∼ σ ∼
1. This case sharesthe same magnetization factor as Case II and III. From the simulation we can see that thegeneral features are similar, but the increases in the poloidal and the toroidal components atthe shock front are nearly identical (Fig. 12 (left), B y and B z ). However, since the emissionregion is a cylinder, a larger fraction (at larger radii) of the emission region possesses theenhanced toroidal component. Hence we expect that the radiation during the flare shouldhave more toroidal contributions. In order to test the effects of the injection rate, we performtwo 3DPol runs: one with an artificially increased injection rate (approximately ten timeshigher than the normal value), and one with no injection; the results are shown in Fig. 12(right). We observe that although the flare level of the high injection case is much higherthan the no injection case, the polarization variations appear similar (Fig. 12 (right)). Thisindicates that the polarization signatures in our calculation are largely due to the magneticfield evolution, instead of the ad hoc nonthermal electron injection rate. We can see that14he PA rotation does not show significant toroidal dominated plateau. A number of flare +PA rotation events feature continuous PA rotations without any obvious plateau step (e.g.,Marscher et al. 2008; Abdo et al. 2010; Chandra et al. 2015). Based on our simulations, wesuggest that such events are likely originating from a sufficiently magnetized environmentwith a relatively fast disturbance traveling through. Also in this case a vigorous restoringphase is present.
5. Discussions and Summary
We have presented the first polarization-dependent radiation modeling based on RMHDsimulations of relativistic shocks in helical magnetic fields in the blazar emission region.We find in Section 3 that the magnetic field evolution during the shock is intrinsicallygoverned by a competition between the shock speed and the magnetization in the emissionregion: the shock/disturbance tries to alter the magnetic topology, while the magnetizationattempts to resist any modifications. We then investigate the parameter space in Section 4 tofurther constrain the shock parameters and the magnetization in the emission environment.In the following we will summarize the major features and illustrate their connections toobservations.We set up an emission region initially pervaded by a purely helical magnetic field withcomparable poloidal and toroidal contributions. Due to the LOS effect, at the beginningthe polarization is dominated by the poloidal component. A flat relativistic disturbance,traveling parallel to B z , will then form a shock wave and propagate through the emissionregion. In a weakly magnetized environment ( σ < σ & p = r c + 2 r c − p is the power-law index and r c is the shock compression ratio mentioned previously.We estimate that the strongest shock in our slow disturbance cases can only reach a non-thermal particle spectrum of power-law index of ∼
4, which is too soft to fit the generalblazar SEDs. In addition, relativistic collisionless shocks in the highly magnetized regimemay be unable to generate fluctuating magnetic fields and hence inhibit Fermi acceleration(e.g., Sironi et al. 2015). Therefore, a slow shock/disturbance may not efficiently acceleratenonthermal particles that are necessary to generate a flare. This poses a tricky problemfor slow shocks. On the other hand, for a relatively fast shock/disturbance, the polariza-tion variations can either fluctuate around some mean value or smoothly rotate in the PAand present a restoring phase at the end of the flare; both are detected in observations(e.g., Larionov et al. 2008; Abdo et al. 2010; Blinov et al. 2015). In summary, we prefer astrongly magnetized emission region, in which fast disturbances can be the driver of theflaring activities.Observations and spectral variability fittings have shown that the typical Lorentz factorof the blazar emission region is within a few tens (e.g., Jorstad et al. 2005; B¨ottcher et al.2013). This is likely to be the upper limit of the shock speed. Therefore, a sufficientlymagnetized emission environment will mostly prohibit any dramatic polarization variations,such as the PA rotations and the sudden changes in the PD (flares in the polarized flux).For the blazar emission regions with a weaker magnetization, the PD flares and the PA rota-tions are triggered by strong shocks. Hence we expect strong nonthermal electron injections.Therefore, intense polarization variations are probably accompanied by vigorous multiwave-length flares. These findings are consistent with observations (Blinov et al. 2015). Anotherimplication is that the blazar emission region presumably maintains a relatively high mag-netization, thus it cannot dissipate too much magnetic energy during the shock-initiatedflaring activities.An alternative may weaken some of the constraints mentioned above. In a highly mag-netized region, magnetic reconnection can also be the driver of flares. Several authors haveillustrated that reconnection can efficiently accelerate nonthermal particles to form a power-law spectrum (Guo et al. 2014; Sironi & Spitkovsky 2014; Li et al. 2015). In our simulations,no preexisting turbulent magnetic field is present; but in reality, some amount of turbulencemay exist. Therefore, during the shock propagation, magnetic reconnection can happen si-multaneously, especially in the compression regions at the shock front. This can provideextra nonthermal particles, thus a slow shock in a partially turbulent, partially ordered16such as a helical geometry) magnetic field structure may generate the necessary nonthermalelectrons. Several models investigating the reconnection in a highly magnetized environmenthave been proposed to explain very fast blazar events (Giannios et al. 2009; Deng et al. 2015;Guo et al. 2015). Meanwhile, recent simultaneous fittings of blazar SEDs, light curves andpolarization signatures also favor a reconnection process (Zhang et al. 2015; Chandra et al.2015). During a reconnection event, the magnetic field topology can be greatly modified,which may give rise to some interesting polarization signatures as well. Therefore, we expectthat even in a highly magnetized environment, if the reconnection drives the flaring activi-ties, considerable polarization variations are possible. Future observations featuring both theradiation and polarization signatures, along with detailed modelings and simulations of theshock and magnetic reconnetion, can possibly distinguish the two mechanisms and furtherconstrain the physics of blazar emission regions.We thank the referee for insightful and constructive comments. HZ thanks Fan Guo,Shengtai Li and Xiaocan Li for helpful discussions. HZ, WD and HL are supported bythe LANL/LDRD program and by DoE/Office of Fusion Energy Science through CMSO.MB acknowledges support by the South African Research Chairs Initiative (SARChI) of theDepartment of Science and Technology and the National Research Foundation of SouthAfrica. Simulations were conducted on LANL’s Institutional Computing machines. REFERENCES
Abdo, A. A., et al., 2010, Nature, 463, 919Achterberg, A., Gallant, Y. A., Kirk, J. G., & Guthmann, A. W. 2001, MNRAS, 328, 393Aharonian, F. A., et al., 2007, ApJ, 664, L71Albert, J., Aliu, E., Anderhub, H., et al. 2007, ApJ, 669, 862Barnacka, A., Moderski, R., Behera, B., Brun, P., & Wagner, S., 2014, A&A, 567, A113Blinov, D., Pavlidou, V., Papadakis, I., et al. 2015, arXiv:1505.07467B¨ottcher, M., & Dermer, C. D., 2010, ApJ, 711, 445 Any opinion, finding and conclusion or recommendation expressed in this material is that of the authorsand the NRF does not accept any liability in this regard.
This preprint was prepared with the AAS L A TEX macros v5.2. xz plane. We only cut part of thesimulation in z direction that is related to the emission region. The emission region rangesfrom − . .
5. The three panels are the toroidal component B y , which represents the B φ in this plane for a helical magnetic field ( B x is trivial in this plane), the poloidal component B z and the velocity of the disturbance V z . All color bars show values in the RMHD codeunits. 21ig. 3.— RMHD simulation results for Case I. The first row shows the magnetic fieldstructures and the plasma velocity when the disturbance is moving into the emission region.The second row shows the same information when the disturbance is moving out. Noticethe color bars for the magnetic fields are different from Fig. 2, and the region is cut smaller.Otherwise plots are as in Fig. 2.Fig. 4.— Hydro parameters in the xz plane for Case I, when the disturbance is movingout of the emission region. The three panels are the plasma density ρ , the pressure P , andvelocity in x direction V x , which represents the velocity in radial. All quantities are in codeunits. 22ig. 5.— RMHD simulation results for Case II. Plots are as in Fig. 3.Fig. 6.— RMHD simulation results for Case III. Plots are as in Fig. 3.23 F l u x =10, =0.01 =10, =1 =3, =1 P D PA Time (d)
Fig. 7.— Time-dependent radiation and polarization signatures. The x axis is time in days.The upper panel shows the relative flux in the optical band. The middle panel is the opticalPD in percentage. The lower panel is the corresponding PA in degree. Black solid lines arefor Case I, red dash for Case II, and blue dot for Case III.24ig. 8.— RMHD simulation results for a weakly magnetized environment ( σ ∼ . V z . The left panel is for a medium speed disturbance(Γ = 6), the right a slow speed disturbance (Γ = 3). Otherwise plots are as in Fig. 3.25 P D Time (d) PA Fig. 9.— Time-dependent radiation and polarization signatures for σ ∼ .
1. The relativeflux is not plotted. Black solid lines are for Γ = 6, red dash are for Γ = 3. Otherwise panelsand curves are as in Fig. 7. 26ig. 10.— RMHD simulation results for a strongly magnetized environment ( σ ∼ P D Time (d) PA Fig. 11.— Time-dependent radiation and polarization signatures for σ ∼
10. Panels andcurves are as in Fig. 9. 28 .00.51.01.52.00102030 0 5 10 15 20 25 3060120180240 F l u x High Injection No Injection P D PA Time (d)
Fig. 12.— RMHD simulations and the associated time-dependent radiation and polarizationsignatures for Γ = 6, σ ∼ ..