Polarization Signatures of Relativistic Magnetohydrodynamic Shocks in the Blazar Emission Region - I. Force-free Helical Magnetic Fields

The blazar emission region is often considered to be an unresolved region near the broad line region of the jet. Generally speaking, the magnetic field topology and evolution as well as plasma 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 emission region 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 can be considered as uncorrelated with the properties of the large-scale jet, even though it is spatially embedded within the jet. On the other hand, the emission region is expected to remain well localized in the jet during flares due to pressure provided by the large-scale jet structure. We assume that flares are due to a disturbance propagating through the emission region. Many observations show that light curves and even the time-dependent polarization signatures appear symmetric in time (e.g., Abdo et al., 2010). This means that the emission region 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 detailed consideration of light-travel-time effects (LTTEs), can naturally explain the apparently time-symmetric light curves and time-dependent polarization signatures. Therefore, we continue to use that assumption here.

Due to the relativistic aberration, even though we are observing blazars nearly along the jet in the observer’s frame (typically, $\theta_{obs,1}\sim 1/\Gamma_{1}$), where $\theta_{obs,1}$ and $\Gamma_{1}$ are the angle between the LOS and the jet direction and the Lorentz factor of the emission region in the observer’s frame, respectively, the angle $\theta_{obs}$ between LOS and the jet axis in the comoving frame is likely around $90^{\circ}$ (if $\theta_{obs,1}\sim 1/\Gamma_{1}$, then $\theta_{obs}\sim 90^{\circ}$). Zhang et al. (2014) have shown that if the LOS is fixed at some other angles, the general trend of polarization variations is only weakly affected. In addition, as is suggested in the bending jet scenario of the polarization signatures, 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 the following simulations, we choose $\theta_{obs}=90^{\circ}$ in the comoving frame, and hence the Doppler factor $\delta\equiv\left(\Gamma_{1}\,[1-\beta_{\Gamma_{1}}\,\cos\theta_{obs,1}]\right)^{-1}\sim\Gamma_{1}$.

We use ideal RMHD simulations to describe the evolution of the system. The underlying model assumes that a plasma jet is traveling with a Lorentz factor of $\Gamma_{0}$ in the observer’s frame. In the comoving frame of the jet, a flow of plasma which is pervaded by a helical magnetic field travels at a Lorentz factor of $\Gamma$. Inside this flow lies a cylindrical emission region containing nonthermal particles. Along the path of the flow it encounters a flat stationary layer of the plasma (the disturbance), forming a shock wave. Thus the total Lorentz factor of the emission region in the observer’s frame is $\Gamma_{1}=\Gamma\times\Gamma_{0}$. Although in the following sections our RMHD simulations use different Lorentz factors $\Gamma$ of the disturbance in the comoving frame of the emission region, $\Gamma_{1}$ of the emission region in the observer’s frame is not constrained. Therefore, we assume a fixed Lorentz factor of $\Gamma_{1}=20$ for the emission 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 is assumed to be in force-balance. The advantage of this magnetic field setup is that it will naturally give rise to comparable poloidal and toroidal contributions, resulting in a relatively low background PD without the need of any turbulence. In addition, this can simplify our setup for the plasma density and the thermal pressure, which are taken to be uniform inside the emission region. In this way, we assume that initially all flow conditions and magnetic fields are laminar, with no preexisting turbulent components. Also, as the RMHD requires an input for the thermal pressure, we assume that in the beginning the plasma is cold and hence the thermal pressure is very small compared to the kinetic energy or the magnetic energy.

In the comoving frame of the emission region, the disturbance and the resulting shock will propagate through the emission region, and change the physical conditions as well as accelerate particles at its location, generating a flare. We apply the simple assumption that the shock will inject fresh nonthermal electrons at the shock front whose energy is equal to a fixed small amount of the local shock kinetic energy. In this way, any deceleration and deformation 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 signatures too 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 the large-scale jet.
4. In the comoving frame of the emission region, the observer is observing from the side of the emission region ($\theta_{\rm 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 comoving frame of the emission region.
9. The disturbance will generate a shock in the emission region. At the shock front, the shock will inject nonthermal particles whose total energy is a constant, small fraction of the local shock kinetic energy.

The above model is realized by the combination of the 3D multi-zone RMHD code LA-COMPASS developed by Li & Li (2003) and the 3D multi-zone polarization-dependent ray-tracing code 3DPol developed by Zhang et al. (2014). Fig. 1 (left) shows a sketch of the simulation setup, and Table 1 presents some major parameters. The RMHD simulation is performed in the comoving frame of the emission region, using Cartesian coordinates. We apply the outflow boundary conditions for the simulation domain; in the following, we will demonstrate that our setup should not suffer from any boundary effects. The simulation domain is pervaded by a helical magnetic field (Fig. 1 (left)). In reality, the poloidal magnetic field lines have to be closed, but based on our assumption the returning magnetic flux is generally outside the simulation domain. The helical magnetic field is assumed to be in force-balance, in the form of

where $J_{0}$ and $J_{1}$ are the Bessel functions of the first kind, $r$ is the radius, and $B_{0}$ and $k$ are normalization 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$ directions and comply with the assumed generally monodirectional $B_{z}$, we add an exponential cutoff for both components when $B_{z}$ approaches zero ($kr\sim 2.4$). However, the exponential cutoff will lead to a nontrivial magnetic force at the edge of the simulation domain, hence we add a plasma pressure to balance it (Fig. 1 (left)).

The emission region is set to be a fixed volume in the simulation domain ($r\sim 2.5$ and $z$ ranges from $-2.5$ to $2.5$, see Fig. 2). Since the simulation domain is in the comoving frame of the emission region, the disturbance will travel at a Lorentz factor of $\Gamma$ in the $z$ direction. We assume that the disturbance is a thin layer ($0.5$ in width) and initially put it at some distance away from the emission region (the disturbance layer ranges from $-8$ to $-7.5$) to allow some time for the shock wave to form in the numerical simulation (see Fig. 1 (left) and Fig. 2). In order to avoid boundary effects in $z$ direction, we use a sufficiently large range in $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,

where $E_{em}=\frac{B^{2}+E^{2}}{8\pi}$ is the electromagnetic energy density, $h=\rho c^{2}+\frac{\hat{\gamma}p}{\hat{\gamma}-1}$ is the specific enthalpy, $\rho$ is the plasma density, $\hat{\gamma}$ is the adiabatic index and $p$ is the thermal pressure. We further assume that the plasma is cold. Hence we choose $\frac{\hat{\gamma}p}{\hat{\gamma}-1}$ to be a small part ($1/30$) of the electromagnetic energy initially. In this way the magnetization factor is approximately $\sigma\sim\frac{E_{em}}{\rho c^{2}}$. The relativistic Alfvén speed can then be expressed as

Since the entire simulation domain is axisymmetric, we present a cut in the $xz$ plane in Fig. 2 to illustrate the initial condition, including the magnetic field and the velocity of the disturbance. 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 $\Gamma=10$; for other $\Gamma$ the velocity 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. It uses a cylindrical geometry, and further divides it evenly into multiple zones in $r$, $\phi$ and $z$ directions. Each zone has its own magnetic field evolution, which is provided by the RMHD simulations. The nonthermal electrons in each zone are assumed to range between a minimal and a maximal Lorentz factor $\gamma_{1}$ and $\gamma_{2}$, with a fixed spectrum of power-law index of $-2$ plus an exponential cutoff at the high-energy end,

The nonthermal electron normalization factor $n$ is assumed to have two components: a fixed uniform initial background density $n_{b}$, and an injected density $n_{i}$ during the shock. As the electron spectrum is fixed, we will only consider the optical light curves and polarization patterns, which do not suffer from any synchrotron-self absorption and Faraday rotation effects for our parameters. As in the RMHD simulation, all calculations are performed in the comoving frame of the emission region. Given the local magnetic field information and the nonthermal electron population, the 3DPol code will calculate the Stokes parameters at various frequencies at every time step in each zone. Next the code will employ ray-tracing to take account of all LTTEs (see Fig. 1 (right) for an illustration), and add up the incoherent Stokes parameters arriving at the observer at the same time, then Lorentz transform to the observer’s frame so as to obtain the total time-dependent radiation and polarization signatures.

We define the PA in the following way. When the electric vector is parallel to the emission region propagation, $PA=0$ (toroidal component dominating). Since the PA has $180^{\circ}$ ambiguity, the toroidal dominance happens at $PA=2N\times 90^{\circ}$, and the poloidal dominance at $PA=(2N+1)\times 90^{\circ}$, where $N$ is an integer.

We first examine a kinetic energy dominated environment, Case I. Such an environment is frequently assumed in a shock-initiated flaring model. We first notice that in a trivially magnetized setting ($\sigma\sim 0.01$), the speed of Alfvén waves, $V_{A}\sim 0.1c$, is much slower than that 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 can expect that the shock/disturbance will strongly compress the plasma at the shock front, and push it along with the propagation. However, in the ideal MHD, the plasma and the field lines are coupled (frozen-in field lines). Therefore, the toroidal magnetic field lines at the shock front 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 to a perfectly ordered poloidal magnetic field (see Fig. 3, $B_{y}$). The physical effects described above are the consequence of a shock in a nearly hydro environment with the assumption of the frozen-in field lines. Therefore, we expect that the above shock modifications to the magnetic field structure is not limited to our initial magnetic field topology, but will happen in 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 pressure and density may be broken. This will lead to some hydrodynamical perturbations in the plasma 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 by Fig. 4 where any plasma velocities in the radial direction that arise from those perturbations are non-relativistic. Therefore, although given enough time they may significantly modify the post-shock structures, they are much slower than the time scale that we are interested in in our RMHD simulation. Therefore, we find that those hydrodynamical effects are of minor importance to our current study, and will not discuss them in the following.

Case II has a moderately magnetized environment ($\sigma\sim 1$). Hence, the magnetic field will actively participate in the shock propagation. Komissarov & Lyutikov (2011) have calculated that the maximal relativistic shock compression ratio in a magnetized environment is given by

where $r_{c}$ and $\sigma^{\prime}$ are the shock compression ratio and the magnetization factor in the shock frame, respectively. As a result, at the shock front, we can see from the simulation that the shock compression becomes weaker due to the higher magnetization (Fig. 5, $B_{y}$). Furthermore, although initially the emission region has a negligible magnetic force, the enhancement in the toroidal component at the shock front will break down this balance, resulting in a contracting magnetic force in the $r$ direction, given by

This contraction is transported by Alfvén waves, which are mildly relativistic in this case. Therefore, in addition to the compression in the $z$ direction at the shock front, there exists a contraction of the plasma in the $r$ direction. Because of the frozen-in field lines, the poloidal component will be strengthened at the shock front (Fig. 5, $B_{z}$). We can see from the simulation that at a later stage, the increase in the poloidal component gradually gets closer to 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 the shock/disturbance, the magnetic field is observed acting to revert to its initial topology (Fig. 5, $B_{y}$ and $B_{z}$).

In the spectral variability fitting, it is frequently argued that during the shock propagation, the local magnetic field either generally maintains its strength and topology (Sokolov et al., 2004; Joshi & Böttcher, 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. Our simulations demonstrate that this can only happen in an adequately magnetized emission region. Therefore, detailed analysis of the magnetic field evolution and shock dynamics should play 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 to some radiation and polarization signatures, but its effect is expected to be weaker than the main shock. Thus in the following we will not consider it.

Now we will consider the radiation and polarization signatures resulting from the above RMHD simulations, calculated using the 3DPol code. The general trends are similar for both cases. Before the shock moves in, the entire emission region is axisymmetric. Although the poloidal and toroidal components are comparable in the force-free setup, due to the LOS effect, the projected toroidal component onto the plane of sky is weaker than the projected poloidal component. Therefore initially the emission region has a total polarization dominated by the poloidal component (Fig. 7). When the shock moves in, the toroidal component is strongly increased and fresh nonthermal electrons are injected at the shock front. However, due to the LTTEs, only a small elliptical region on the right that is near the observer is seen (the near side, see Fig. 1 (right), red). Therefore, the flux is seen to gradually increase. Nevertheless, since the toroidal enhancement is strong, a small flaring region is adequate to push the polarization to be dominated by the toroidal component. As a result, the PA quickly rotates to $180^{\circ}$ and forms a plateau, which corresponds to the toroidal domination. The PD first experiences a drop which indicates a switch in the domination between the two components, then climbs up to a high level, implying a strong toroidal dominance (Fig. 7). When the shock is about to leave the emission region, the flaring region reaches maximum (Fig. 1 (right), green). Hence the flux peaks. After that, the flaring region moves far from the observer to the left (the far side, Fig. 1 (right), blue), so the flux gradually decreases in an apparently symmetric pattern. However, as some plasma and the coupled magnetic field lines have been pushed away by the shock/disturbance, the total magnetic field strength is smaller than the initial state. Hence the flux level is lower than the initial value. For the same reason, the toroidal contribution is weakened near the end of the flare, hence the PD rises up higher than the initial value, revealing a stronger poloidal contribution. Since the helical magnetic field direction on the far side is opposite to that on the near side, the PA instead completes a $180^{\circ}$ rotation to the initial value ($180^{\circ}$ PA ambiguity, Fig. 7).

The major differences between Case I and II are the following. First, due to the stronger magnetic field, Case II has a weaker compression in the toroidal component at the shock front than Case I, while it experiences a rise in the poloidal contribution. Therefore, the PA rotation in Case II appears smoother and has a shorter toroidal dominated plateau than Case I. Also the maximal flare level in Case II is lower. More importantly, near the end of the flare, the moderate magnetization in Case II is acting to restore the initial magnetic topology, 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 shown that at the end of a flare + PA swing event in 3C 279, the PD rises up above the initial value, which is immediately followed by a significant restoration to approximately the initial PD level. The same restoring phases are seen in a number of flare + PA rotation events as well (Larionov et al., 2008; Larionov et al., 2013; Morozova et al., 2014). On the other hand, in Case I the magnetization is too weak, thus the PD rises up to $\sim 70\%$ at the end of the flare and shows no restoration, indicating a purely poloidal magnetic structure. As is mentioned above, such a phenomenon is generally the consequence of a trivially magnetized environment, and it is expected to happen with other initial magnetic field topologies. Therefore, we argue that such 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 by examining Case III, which shares the same magnetization as Case II but has a slower disturbance. A slower disturbance has less kinetic energy. Therefore, as the shock/disturbance propagates, it cannot exert a very strong pressure onto the plasma; in return the plasma will decelerate the shock/disturbance. Thus we see in Fig. 6 that the shock/disturbance is considerably slower than the initial value ($\Gamma=3$). Therefore, the flare duration is much longer than the previous cases (Fig. 7). Also the shock/disturbance will provide less compression in the plasma at the shock front. As a result, we observe in the simulation that the toroidal component is only moderately increased and stretched at the shock front and downstream of the shock/disturbance, respectively (Fig. 6, $B_{y}$). Moreover, the shower shock/disturbance takes longer time to propagate through the emission region, allowing for further information exchange by Alfvén waves. Therefore, the contraction of the plasma at the shock front, due to the enhanced toroidal component, is seen to catch up with the shock compression. This leads 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 contribution. Therefore, the PD still displays a decrease at the middle of the flare (Fig. 7). Meanwhile, downstream of the shock/disturbance, Alfvén waves restore the initial magnetic field topology to a large extent (Fig. 6, $B_{y}$ and $B_{z}$). Hence the final PD is not too much higher than the initial value, and gradually recovers after the flare (Fig. 7). Another interesting effect is that due to the weaker kinetic energy in the shock/disturbance and longer influence of Alfvén waves, the shape of the shock/disturbance is distorted (Fig. 6, $V_{z}$). Thus the polarization patterns, especially the PA, appear less symmetric in time than the previous cases (Fig. 7).

We study a weakly magnetized setup ($\sigma\sim 0.1$), with the application of moderate ($\Gamma=6$) and slow ($\Gamma=3$) disturbances. Due to the small magnetization factor, the mediation by Alfvén waves is slow. Nevertheless, especially in the case of a slow disturbance, the effects of the magnetization cannot be overlooked: the poloidal component is slightly enhanced at the shock front, and the shape of the shock/disturbance is changed (Fig. 8, $B_{z}$ and $V_{z}$). As a consequence, the polarization patterns appear relatively smooth, in particular the PA rotation (Fig. 9). However, downstream of the shock/disturbance, Alfvén waves fail to restore the initial magnetic topology. Hence at the end of the flare, the PD rises up and maintains a high level ($\sim 40\%$), without showing any restoring phase (Fig. 9). This implies that in a weakly magnetized emission environment, shocks can easily alter the magnetic topology, causing significant changes in the polarization signatures. On the contrary, such polarization variations are seldom reported in observations. Therefore, we do not favor a weakly magnetized emission environment.

Next we consider an emission region with a high magnetization factor ($\sigma\sim 10$). Here Alfvén waves become relativistic, thus the shock compression is much weaker. All effects of the magnetization as in Case II and III show up. The difference is that even in the case of a medium speed disturbance ($\Gamma=6$), the increase in the poloidal component is higher than that in the toroidal field at the shock front (Fig. 10, $B_{y}$ and $B_{z}$). Therefore, the PA patterns for both cases exhibit no rotation (Fig. 11). Additionally, the shock/disturbance is greatly distorted. Hence both the PD and the PA patterns appear rather erratic, showing a lot of bumps throughout the flare (Fig. 11), even though the general geometry is rather symmetric. 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 recovers at the end of the flare (Fig. 11).

Notice that for a sufficiently fast shock/disturbance, a PA rotation may reappear. However, we expect that such a rotation is likely to be very smooth, without a long toroidal dominated plateau as in Case I to III. Meanwhile the PD variation can be within some small value. This may agree with the observed flare + PA rotation events. To summarize, a strongly magnetized emission environment may generate the typical erratic perturbations in the polarization signatures in the condition of a relatively slow disturbance, together with the smooth flare + PA rotation events provided a relatively fast disturbance. Thus we prefer an emission region with high magnetization.

Finally we take a look at an intermediate condition, $\Gamma=6$ and $\sigma\sim 1$. This case shares the same magnetization factor as Case II and III. From the simulation we can see that the general features are similar, but the increases in the poloidal and the toroidal components at the shock front are nearly identical (Fig. 12 (left), $B_{y}$ and $B_{z}$). However, since the emission region is a cylinder, a larger fraction (at larger radii) of the emission region possesses the enhanced toroidal component. Hence we expect that the radiation during the flare should have more toroidal contributions. In order to test the effects of the injection rate, we perform two 3DPol runs: one with an artificially increased injection rate (approximately ten times higher 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 higher than the no injection case, the polarization variations appear similar (Fig. 12 (right)). This indicates that the polarization signatures in our calculation are largely due to the magnetic field evolution, instead of the ad hoc nonthermal electron injection rate. We can see that the 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, we suggest that such events are likely originating from a sufficiently magnetized environment with a relatively fast disturbance traveling through. Also in this case a vigorous restoring phase is present.