© The Authors 2025

1. Introduction

Multifrequency polarimetry is emerging as a powerful probe of blazar jets, especially with the advent of the Imaging X-ray Polarimetry Explorer (IXPE) satellite (Weisskopf et al. 2016). IXPE has opened a new window for polarimetry in the X-ray band, providing unprecedented insights into the behavior of these extreme astrophysical objects. Blazars at the low-power end of the blazar sequence, i.e., the most efficient particle accelerators, are excellent targets for IXPE because their emission peaks in the X-ray band (Tavecchio 2021). This is particularly true for High-frequency-peaked BL Lacs (HBLs) and Extreme High-frequency-peaked BL Lacs (EHBLs), where both optical and X-ray emissions can be attributed to synchrotron radiation from a population of nonthermal electrons.

Several multiwavelength polarimetric campaigns have been organized to measure the polarization properties of HBLs and EHBLs across the electromagnetic spectrum (e.g., Liodakis et al. 2022; Di Gesu et al. 2022, and Ehlert et al. 2023). Observations of blazars in a quiescent state, where the polarization degree and the electric vector position angle (EVPA) do not show significant temporal variability, have revealed the following key characteristics. First, the polarization degree increases with frequency, a phenomenon known as chromaticity. Specifically, the observed X-ray polarization degree is typically Π_x∼10−20%, which is significantly higher than the optical and radio polarization degree, with Π_x/Π_o∼2−6, where Π_o is the optical polarization degree. Such strong chromaticity suggests that the mechanisms influencing polarization are energy-dependent. Second, the EVPA is weakly chromatic, with the X-ray and optical EVPAs being approximately aligned (ψ_x∼ψ_o), indicating a stable magnetic field orientation in the emitting region. Moreover, the EVPA tends to align with the projection of the jet axis on the plane of the sky. However, this alignment is not consistent across all sources and comes with certain caveats (more details in the discussion). Finally, IXPE also observed an episode of temporal variability, specifically a rotation of the X-ray polarization angle observed in the HBL Mrk421 in June 2022 by more than 360° in 5 days (Di Gesu et al. 2023).

Several models have been proposed to explain the recent results on the multiwavelength polarization of HBLs and EHBLs. For instance, some models are based on the global structure of the jet magnetic field (e.g., Bolis et al. 2024). Currently, the most accepted scenario is the energy-stratified shock model (Angelakis et al. 2016; Tavecchio et al. 2018). Supposing that nonthermal electrons are accelerated by a shock, those emitting X-ray photons cool down in a small region close to the shock front. In this region, the magnetic field is ordered by shock compression and there is a strong self-generated component, implying a large degree of polarization. On the other hand, the particles emitting optical and radio photons cool down in a larger region and are advected by the flow farther downstream. In this case, electrons are subject to several radiation zones with different magnetic field orientations, implying a low polarization degree. A likely origin for the shock responsible for blazar emission is the recollimation of the jet by external gas (Tagliaferri et al. 2008; Marscher 2014; Zech & Lemoine 2021; Tavecchio et al. 2022; Costa et al. 2024) and this is the scenario that we adopt in the following.

Our goal is to give a more quantitative description of the energy-stratified model in the framework of a recollimating jet scenario, starting directly from relativistic magnetohydrodynamic (RMHD) simulations. From fluid simulations, we obtained the large-scale structure of the magnetic and velocity field, necessary to estimate the nonthermal particle emission and its beaming. However, fluid simulations cannot resolve kinetic scales and provide no information about the injection and acceleration of nonthermal particles, which requires tracking particles on scales of the electron gyroradius, as done in particle-in-cell (PIC) simulations.

One possible solution to bridge the gap between fluid and kinetic scales is the Lagrangian macroparticle approach. A macroparticle is an ensemble of nonthermal particles decoupled from the fluid. The macroparticles move along the streamlines, and their spectra are updated using the local fluid conditions. The particles sample the fluid, and their integrated emission is used to estimate the radiative output of the jet (Mimica et al. 2009; Mimica & Aloy 2012). Relativistic hydrodynamic simulations have frequently been considered because of their lower computational cost. In this case, the magnetic field is obtained assuming an equipartition between the magnetic and internal energy density. Moreover, it was standard to consider a fixed power law index for the spectra of the macroparticle injected at the shock, neglecting the dependence of the acceleration process on the shock features.

These two limitations have been surpassed by the approach outlined in Vaidya et al. (2018). This hybrid framework estimates the emission directly from RMHD simulations and implements a sub-grid model for shock acceleration. The magnetic field is used to calculate the synchrotron losses and emission of the macroparticles, accounting for beaming and the relativistic aberration of synchrotron polarization. Moreover, when macroparticles cross a shock, their spectra are updated using the local fluid conditions, specifically the magnetic orientation and the shock compression ratio.

In this paper, we implemented a new post-processing code, based on the work of Vaidya et al. (2018) and Mukherjee et al. (2021), which can calculate the radiative output of stationary RMHD simulations. Comparing our results with early PIC simulations, we find that shocks formed by jet recollimation are predominantly superluminal, which restricts particle acceleration in a laminar flow (Crumley et al. 2019 for the transrelativistic case, Sironi & Spitkovsky 2009, 2011 for the relativistic case). However, recent PIC simulations indicate that acceleration can still occur in a superluminal configuration if small-scale turbulence or inhomogeneities are present (Bresci et al. 2023; Demidem et al. 2023). With particle acceleration turned on, we focused on investigating the energy-stratified shock scenario, specifically evaluating whether it can replicate the polarimetric features for HBLs and EHBLs detected by IXPE and the simultaneous multiwavelength campaigns. In Sect. 2 we describe the RMHD simulations and in Sect. 3 the post-processing code. The application is reported in Sect. 4 and in Sect. 5 we discuss the results.

2. RMHD setup

2.1. Context

We simulated the recollimation of a relativistic jet at parsec-scale distance from the supermassive black hole, where the high-energy emission region is probably located (Tavecchio et al. 2010a). Aside from occasional flares, the observed emission remains stable throughout day-long observation periods. This level of stability would be difficult to explain if the emission were generated by a shock traveling through the jet. Assuming that recollimation and reflection shocks are responsible for particle acceleration is not a standard choice since in blazar modeling it is often postulated that the emission region is downstream of a shock of normal incidence traveling down the jet. Possible physical realizations of this setup could be represented by a dense blob of plasma overtaking the jet plasma or by internal shocks developing between portions of the flow characterized by different speeds (Rees 1978; Spada et al. 2001). However, a difficulty of such a setup (already pointed out by e.g., Tagliaferri et al. 2008) is that in an observed time interval Δt, the shock moves by a distance $\mathrm{\Delta }z\sim c\mathrm{\Delta }t{\mathrm{\Gamma }}_{\mathrm{sh},\mathrm{obs}}^2$, where Γ_sh,obs is the Lorentz factor of shock in the observer frame. Assuming typical values, namely Δt = 1 day and Γ_sh,obs = 30, the distance amounts to ∼4 pc. During the propagation (and the consequent expansion) the physical parameters of the plasma in the jet (density, magnetic field) would drastically change, greatly modifying the emission properties (Blandford & Königl 1979). Such a conclusion is in patent disagreement with the stability of the emission during the considered observation periods (e.g., Liodakis et al. 2022). We remark that the large Lorentz factor of the shock in the observer frame is a direct consequence of the fact that to properly model blazar emission we need to Doppler boost the emission from the downstream plasma, while for a relativistic shock of normal incidence, the velocity of the downstream flow in the shock frame is necessarily sub-relativistic, hence the shock frame must move relativistically in the observer frame. This problem can be overcome in the case of an oblique shock, since in that case only the normal component of the downstream velocity must be sub-relativistic with respect to the shock. In this case, the shocked plasma can maintain a substantial bulk Lorentz factor with respect to the observer even if the shock is stationary in that frame. The low temporal variability can be related to some permanent and fixed (in the observer frame) structure in the jet, an oblique standing recollimation region, shaped by the pressure contrast between the jet and the external medium, e.g., Komissarov & Falle (1997) and Bodo & Tavecchio (2018). Moreover, recollimation shocks could also be a natural outcome of MHD jets even without external confinement, e.g., Jannaud et al. (2023).

We started from the axisymmetric hydrodynamic setup described in Costa et al. (2024), but we included the magnetic field. We considered small magnetization values, for which the jet is expected to be unstable (Matsumoto et al. 2021; Boula et al., in prep.) against recollimation instabilities. However, the first series of recollimation-reflection shocks is not disrupted (Costa et al. 2024).

2.2. Numerical framework

Following Costa et al. (2024), we looked for an axisymmetric steady-state solution by performing relativistic magnetohydrodynamic simulations until stationarity was reached. We used the “RMHD” module of the state-of-the-art PLUTO code for astrophysics (Mignone et al. 2007, 2012), which solves the equations of ideal RMHD

$$\frac{\partial }{\partial t}\left(\begin{array}{cc}d\hfill & \hfill \\ \text{m}\hfill & \hfill \\ e_t\hfill & \hfill \\ \text{B}\hfill & \hfill \end{array}\right)+\nabla ·\left(\begin{array}{cc}d\text{v}\hfill & \hfill \\ {\omega }_t{\mathrm{\Gamma }}^2\text{v}\text{v}-\text{b}\text{b}+p_t\text{I}\hfill & \hfill \\ \text{m}\hfill & \hfill \\ \text{v}\text{B}-\text{B}\text{v}\hfill & \hfill \end{array}\right)=\left(\begin{array}{cc}0\hfill & \hfill \\ {\text{f}}_g\hfill & \hfill \\ {\text{f}}_g·\text{v}\hfill & \hfill \\ 0\hfill & \hfill \end{array}\right),$$(1)

where the speed of light c is set to 1, and a factor $\sqrt{4\pi }$ is reabsorbed in the definition of B. The set of primitive variables consists of the rest frame density, the thermal pressure, the three-velocity in the lab frame, and the magnetic field three-vector in the lab frame (ρ,p,v,B). Γ is the Lorentz factor, f_g is the force density in the lab-frame, and I is the unit 3×3 tensor. The conserved quantities are B, and the laboratory frame mass, momentum, and energy densities, respectively defined as

$$\begin{array}{c}\hfill d=\rho \mathrm{\Gamma }\hfill \\ \hfill \text{m}=w_t{\mathrm{\Gamma }}^2\text{v}-b^0\text{b}\hfill \\ \hfill e_t=w_t{\mathrm{\Gamma }}^2-b^0{b}^0-p_t\hfill \end{array},\text{with}\begin{array}{c}\hfill b^0=\mathrm{\Gamma }\text{v}·\text{B}\hfill \\ \hfill \text{b}=\text{B}/\mathrm{\Gamma }+\mathrm{\Gamma }\left(\text{v}·\text{B}\right)\text{v}\hfill \\ \hfill w_t=\rho h+{\text{B}}^2/{\mathrm{\Gamma }}^2+{\left(\text{v}·\text{B}\right)}^2\hfill \\ \hfill p_t=p+\frac{{\text{B}}^2/{\mathrm{\Gamma }}^2+{\left(\text{v}·\text{B}\right)}^2}{2}\hfill \end{array}.$$(2)

The system of RMHD equations is closed by the Taub-Matthews equation of state, defining the rest frame relativistic specific enthalpy as (Mignone et al. 2005)

$$h=\frac{5}{2}\mathrm{\Theta }+\sqrt{\frac{9}{4}{\mathrm{\Theta }}^2+1},$$(3)

with Θ=p/(ρc²), approximating the Synge EoS of a single-specie relativistic perfect fluid (Synge & Morse 1958).

The numerical methods included a linear reconstruction scheme (Mignone et al. 2007), second-order Runge-Kutta time integration (Butcher 2000), and the HLLD Riemann solver (Mignone et al. 2009). We adopted a constrained transport approach, which ensures that the divergence of the magnetic field is null throughout the simulation (Mignone & Del Zanna 2021).

The simulation was carried out in 2D cylindrical coordinates (R,z)∈[0,6]×[1,30] in units of z₀, defined below, with 1400×2200 points. Boundary conditions were set to outflow at all boundaries except at the axis, where they were reflective, and at the base, at z=z₀. There, analytical profiles were defined constant to continuously inject the jet at radii R<R_j=θ_jz₀, and fixed the environment profile at larger radii too. The steady state obtained from the final output, at t = 3000t₀, where t₀=z₀/c, was finally interpolated to a 3D Cartesian domain (x,y,z)∈[−0.5,0.5]×[−0.5,0.5]×[1,10.5] with 500×500×750 points uniformly spaced.

2.3. Parameters and magnetic field

We simulated a recollimating jet, following the setup described in Costa et al. (2024). The hydrodynamic setup was the same: an expanding relativistic conical jet of width θ_j was initialized in the computational domain at a distance z₀ from the cone vertex, surrounded by a confining environment. In these simulations, we further included the magnetic field.

The magnetic field is helical, defined in cylindrical coordinates in the laboratory frame with components

$$B_R=\frac{\alpha B_0}{r^2}{\text{e}}^{-{\mathcal{X}}^2}\frac{R}{r},$$(4)

$$B_z=\frac{\alpha B_0}{r^2}{\text{e}}^{-{\mathcal{X}}^2}\frac{z}{r},$$(5)

$$B_{\varphi }=\mathrm{\Gamma }\frac{B_0}{r}\sqrt{{\text{e}}^{-2{\mathcal{X}}^2}-\frac{{\psi }_{\chi }}{2{sin}^2\theta }\left[{\psi }_{\chi }-{\psi }_{\chi }{\text{e}}^{-2{\mathcal{X}}^2}+\sqrt{2\pi }\text{erf}\left(\sqrt{2}\mathcal{X}\right)\right]},$$(6)

where $r=\sqrt{R^2+z^2}$, and $\mathcal{X}=\left(cos\theta -1\right)/{\psi }_{\chi }$. Γ and θ=cos⁻¹(z/r) are functions of the position, α and B₀ are parameters related to the pitch and to the amplitude of the field in the rest frame; finally, ψ_χ is the scale length for the decay of the poloidal field, which is assumed to have a Gaussian profile. In our simulations, we set ψ_χ = 10⁻²θ_j, to make the field decay to 0 outside the jet.

Inside the computational domain, the poloidal field $B_{\mathrm{pol}}=\sqrt{B_R^2+B_z^2}$ peaks in (0,z₀) to αB₀, while the toroidal field B_tor=B_ϕ peaks approximately in (0.5θ_j,z₀) to max(B_tor) = 0.7ΓB₀.

These profiles ensure that the magnetic field is in an equilibrium configuration at the lower boundary when α = 1. In case α<1, the jet thermal pressure is defined to compensate for the lower poloidal field,

$$p_j=p_{j,\mathit{HD}}+(1-{\alpha }^2)\frac{B_0^2}{2}\frac{{\text{e}}^{-2{\mathcal{X}}^2}}{r^{2{\gamma }_{\text{ad}}}},$$(7)

where p_j,HD is the hydrodynamic adiabatic pressure profile used in Costa et al. (2024) and γ_ad = 5/3 is the adiabatic index. More details on the magnetic field prescription can be found in the Appendix A.

To simulate the properties of HBLs, we constrained the emitting region to be under/around the parsec scale, where the bulk Lorentz factor is expected to be on the order of 10, the proper frame magnetic field (or better, the component responsible for the emission) around 10⁻² G, and the power around 10⁴³ erg/s (e.g., Tavecchio et al. 2010b). The jet steady state is determined by the environment, and by its own properties, parameterized by a few quantities that we defined at the injection boundary, which we set at z₀ = 0.1 pc. In general, we use the suffixes “j” and “ext” respectively to refer to jet and environment parameters. We defined Γ_j = 10, the jet's opening angle θ_j = 0.1, the jet-enviornment density contrast ν=ρ_j(z₀)/ρ_ext(z₀) = 7.6×10⁻⁶, and we considered a cold jet, with h_j(z₀)≃1 (the jet-environment pressure ratio is p_r=p_j(z₀)/p_ext(z₀) = 10⁻³).

The confining external medium is at rest in the laboratory frame, it is isothermal, with T = 1.5×10⁷ K, and it is stratified, following a power law pressure profile, $p_{\text{ext}}(z)=p_{\text{ext}}(z_0){\left(z/z_0\right)}^{-0.5}$, where p_ext(z₀) = 3×10⁻⁶ρ_ext(z₀)c² and ρ_ext(z₀) = 10⁵ m_p cm⁻³, was chosen to tune the jet power to about 10⁴³ erg/s. The environment is kept in hydrodynamic equilibrium with a force density f_g=∇p_ext.

For the magnetic field, we set the rest frame jet magnetization at the nozzle, defined as

$${\sigma }_j=\frac{B_0^2}{{\rho }_j(z_0)c^2}={10}^{-2},$$(8)

so the proper-frame amplitude of the magnetic field at jet injection is

$$B_0=1.2\times {10}^{-2}\left(\frac{{10}^{-2}}{\sigma }\right)\left(\frac{{\rho }_{\text{ext}}(z_0)}{{10}^5{m}_p{\text{cm}}^{-3}}\right)\text{G}.$$(9)

We studied three cases at different pitches: Case A assumes a force-free field, with α = 1, Case B considers a toroidal-dominated field, with α = 0.5, and Case C, where α = 0, is for a purely toroidal field.

We finally observed that, since σ_j≪1 in all cases we considered, the Poynting flux contribution to the jet power L_j was negligible with respect to the hydrodynamic jet power, L_j,HD, and we found

$$L_j\simeq L_{j,\mathit{HD}}=(\pi z_0^2{\theta }_j^2){\rho }_j{h}_j{c}^2{\mathrm{\Gamma }}_j^2{v}_j$$(10)

$$\simeq {10}^{43}{\left(\frac{{\theta }_j}{0.1}\right)}^2\left(\frac{\nu }{{10}^{-5}}\right)\left(\frac{{\rho }_{\text{ext}}(z_0)}{{10}^5{m}_p{\text{cm}}^{-3}}\right)\left(\frac{h_j}{1}\right){\left(\frac{{\mathrm{\Gamma }}_j}{10}\right)}^2\text{erg s}{}^{-1}.$$(11)

3. Post-processing calculation of particle evolution and emission

As mentioned above, the emission properties of the recollimation shock structure were evaluated in post-processing, tracking particle evolution and emission by means of a Lagrangian macroparticles approach.

A Lagrangian macroparticle is a cloud of actual particles (electrons, protons, etc.) represented by a Dirac delta in physical space and characterized by a nontrivial distribution in phase space. The Lagrangian macroparticles are advected in the Eulerian grid according to

$$\frac{d{x}_p}{\mathit{dt}}=v_p(x_p),$$(12)

where v_p(x_p) is the fluid velocity interpolated at the macroparticle position x_p. In our code, the fluid quantities were interpolated at macroparticle positions using linear interpolation. Eq. (12) was numerically solved by the second-order Runge-Kutta scheme.

While a macroparticle is moving in the fluid, its spectrum is simultaneously updated by using the local fluid conditions. The spectral time evolution is governed by the relativistic transport equation. For the full equation under the assumption of momentum isotropy, refer to Webb (1989). For our application, we can neglect spatial diffusion, shear acceleration, second-order Fermi acceleration, and non-inertial energy change. The transport equation reduces to

$$p^2\frac{\mathit{df}}{d\tau }+\frac{\partial }{\partial p}\left[-\frac{p^3}{3}f{\nabla }_{\mu }{u}^{\mu }+\langle \dot{p}{\rangle }_{\text{rad}}f\right]=-p^{2}f{\nabla }_{\mu }{u}^{\mu },$$(13)

where u^μ is the fluid four-velocity, while all other quantities are expressed in the fluid rest frame. We denote the isotropic phase-space distribution function by f and the Lagrangian derivative with respect to proper time by d/dτ=u^μ∇_μ. The two terms in the square brackets respectively describe the adiabatic and radiative losses ($\langle \dot{p}{\rangle }_{\text{rad}}$ is the average momentum loss due to radiative mechanisms). The right-hand side term accounts for changes in macroparticle number density due to fluid expansion or contraction. When considering a single macroparticle, the total number of particles remains constant. However, as the surrounding fluid expands or contracts, its volume changes, leading to a variation in number density. We define $\mathcal{N}(p,\tau )=\int d\mathrm{\Omega }{p}^{2}f=4\pi p^{2}f$, the number of particles per unit of volume and momentum. Since nonthermal particles are ultrarelativistic, their energy can be expressed as E=pc, which in turn implies $\mathcal{N}(E,\tau )=\mathcal{N}(p,\tau )/c$, where $\mathcal{N}(E,\tau )$ is the number of particles per unit volume and energy. Integrating Eq. (13) over the solid angle and moving to the energy space gives

$$\frac{d\mathcal{N}(E,\tau )}{d\tau }+\frac{\partial }{\partial E}\left[\left(-\frac{E}{3}{\nabla }_{\mu }{u}^{\mu }+\dot{E}\right)\mathcal{N}(E,\tau )\right]=-\mathcal{N}(E,\tau ){\nabla }_{\mu }{u}^{\mu },$$(14)

where ${\dot{E}}_r=\langle \dot{p}{\rangle }_{\text{rad}}/p^2$. Finally, it is possible to eliminate the term on the right side by introducing $\chi =\mathcal{N}(E,\tau )/n$, where n is the fluid number density. Using particle number conservation, Eq. (14) can be rewritten as

$$\frac{d\chi }{d\tau }+\frac{\partial }{\partial E}\left[\left(-\frac{E}{3}{\nabla }_{\mu }{u}^{\mu }+{\dot{E}}_r\right)\chi \right]=0.$$(15)

The first term in the round brackets describes the energy change due to adiabatic expansion or compression, while the second term accounts for radiative losses, which, at the moment, includes synchrotron losses exclusively,

$${\dot{E}}_r=\frac{4{\sigma }_{T}c{\beta }_e^2{U}_B}{3m_e^2{c}^4}{E}^2,$$(16)

where σ_T is the Thomson cross section, β_e≈1 is the electron velocity in units of c, m_e is the electron mass, and U_B=B²/8π is the fluid rest frame magnetic energy density (in contrast with the previous Section, we do not reabsorb any factor $\sqrt{4\pi }$ in the field definition from now on). For simplicity, we neglect inverse-Compton losses, a suitable choice for sources like HBLs.

Eq. (15) is a first-order hyperbolic partial differential equation with variable coefficients and is numerically solved in our code using a second-order implicit-explicit Runge-Kutta scheme, specifically the strong stability preserving algorithm (2,2,2), described in Pareschi & Russo (2005) (see Appendix B for more details). For self-consistency, the nonthermal particles associated with a macroparticle must remain within the computational cell that contains the macroparticle itself. This requires the maximum Larmor radius of these nonthermal particles to be smaller than the size of the computational cell. Consequently, we imposed an upper limit on the energy of the nonthermal particles (Vaidya et al. 2018)

$$E_{\max }=\frac{\mathit{eB}min(\mathrm{\Delta }x,\mathrm{\Delta }y,\mathrm{\Delta }z)}{2{\beta }_e},$$(17)

where e is the electron charge, B is the laboratory-frame magnetic field interpolated at the macroparticle position, while Δx, Δy, and Δz are the cell dimensions. This energy limit was also enforced during the shock update procedure.

The macroparticle spectra are necessary to calculate their radiative output, which corresponds solely to synchrotron emission for our purposes. For easier comprehension, in the following the quantities marked with a prime are described in the fluid rest frame, to distinguish them from the equivalent quantities defined in the laboratory frame. The synchrotron specific emissivity of a single macroparticle is

$$J_{\mathrm{syn}}^{\prime }({\nu }^{\prime })=\frac{\sqrt{3}{e}^3}{4\pi m_e{c}^2}|B^{\prime }\times {\hat{\text{n}}}^{\prime }|\int {\mathcal{N}}^{\prime }(E^{\prime })F\left(x\right){\mathit{dE}}^{\prime }.$$(18)

where ν′ is the fluid rest frame frequency and ${\hat{\text{n}}}^{\prime }$ is the versor parallel to the line of sight. Similarly, the linearly polarized specific emissivity of a single macroparticle is given by

$$J_{\mathrm{pol}}^{\prime }({\nu }^{\prime })=\frac{\sqrt{3}{e}^3}{4\pi m_e{c}^2}|B^{\prime }\times {\hat{\text{n}}}^{\prime }|\int {\mathcal{N}}^{\prime }(E^{\prime })G\left(x\right){\mathit{dE}}^{\prime }.$$(19)

The two functions in the integral are defined as

$$F(x)=x{\int }_x^{\infty }{K}_{5/3}(t)\mathit{dt},$$(20)

$$G(x)=x{K}_{2/3}(x),$$(21)

$$\text{with}\hspace{0.5em}x=\frac{4\pi m_e^3{c}^5{\nu }^{\prime }}{3e{E}^{\prime }|B^{\prime }\times {\hat{\text{n}}}^{\prime }|},$$(22)

where K_a is the modified Bessel function of order a (Ginzburg & Syrovatskii 1965). The emissivities are expressed in the fluid rest frame, as the other quantities. However, we are interested in the emissivity in the observer frame and the RMHD simulations exclusively provide us with quantities in the laboratory frame. Therefore, first we have to transform the magnetic field and line-of-sight direction from the laboratory to the fluid rest frame,

$${\hat{\text{n}}}^{\prime }=D\left[\hat{\text{n}}+\left(\frac{{\mathrm{\Gamma }}^2}{\mathrm{\Gamma }+1}\beta ·\hat{\text{n}}-\mathrm{\Gamma }\right)\beta \right],$$(23)

$$B^{\prime }=\frac{1}{\mathrm{\Gamma }}\left[B+\frac{{\mathrm{\Gamma }}^2}{\mathrm{\Gamma }+1}(\beta ·B)\beta \right],$$(24)

where Γ is the fluid bulk Lorentz factor and $D=[\mathrm{\Gamma }(1-\beta ·\hat{\text{n}}){]}^{-1}$ is the relativistic Doppler factor. The emissivities are then calculated in the fluid rest frame and finally boosted to the observer frame, whose quantities are indicated with a tilde,

$${\tilde{J}}_{\text{syn}}(\nu )=D^2{J}_{\text{syn}}^{\prime }({\nu }^{\prime }),$$(25)

$${\tilde{J}}_{\text{pol}}(\nu )=D^2{J}_{\text{pol}}^{\prime }({\nu }^{\prime }),$$(26)

where $\tilde{\nu }=D{\nu }^{\prime }$ is the frequency in the observer frame. This procedure is repeated for each macroparticle, and their emissivities are deposited on the fluid simulation grid using the second-order accurate triangular-shaped cloud method, a standard technique also employed in PIC simulations (Birdsall & Langdon 1991). Next, we introduce a new grid, with the Z axis along the line of sight, while the X and Y axes on the sky plane. After interpolating the emissivities from the simulation to the new grid, the synchrotron specific intensity I and flux F are computed,

$$\tilde{I}(\nu ,X,Y,\theta )=\int {\tilde{J}}_{\mathrm{syn}}(\nu ,X,Y,Z)\mathit{dZ},$$(27)

$$\tilde{F}(\nu ,\theta )=\iint \tilde{I}(\nu ,X,Y)\frac{\mathit{dXdY}}{{\mathcal{D}}^2},$$(28)

where $\mathcal{D}$ is the source distance and θ is the observer viewing angle. In the end, we calculate the synchrotron polarization degree properties. For that, we need the Stokes parameters Q and U (circular polarization is neglected), but relativistic effects must be taken into account. Following Lyutikov et al. (2003) and Del Zanna et al. (2006), the two Stokes parameters on the sky plane are given by

$$\tilde{Q}(\nu ,X,Y,\theta )=\int {\tilde{J}}_{\text{pol}}(\nu ,X,Y,Z)cos2\psi \mathit{dZ}$$(29)

$$\tilde{U}(\nu ,X,Y,\theta )=\int {\tilde{J}}_{\text{pol}}(\nu ,X,Y,Z)sin2\psi \mathit{dZ}$$(30)

where ψ is the EVPA or polarization angle, and its sine and cosine for a specific cell are equal to

$$cos2\psi =\frac{q_X^2-q_Y^2}{q_X^2+q_Y^2},sin2\psi =-\frac{2q_X{q}_Y}{q_X^2+q_Y^2},\text{with}\hspace{0.5em}$$(31)

$$q_X=(1-{\beta }_Z)B_X-{\beta }_X{B}_Z,q_Y=(1-{\beta }_Z)B_Y-{\beta }_Y{B}_Z,$$(32)

where β_X,Y,Z and B_X,Y,Z are the interpolation of the velocity and magnetic field on the new grid. The polarization degree Π and angle (measured clockwise from the Y axis) on the sky plane as functions of the frequency are given by

$$\tilde{\mathrm{\Pi }}(\nu ,X,Y,\theta )=\frac{\sqrt{\tilde{Q}(\nu ,X,Y,\theta {)}^2+\tilde{U}(\nu ,X,Y,\theta {)}^2}}{\tilde{I}(\nu ,X,Y,\theta )},$$(33)

$$tan2\tilde{\psi }(\nu ,X,Y,\theta )=\frac{\tilde{U}(\nu ,X,Y,\theta )}{\tilde{Q}(\nu ,X,Y,\theta )}·$$(34)

Since the source location is unresolved, we integrate the Stokes parameters over the sky plane, obtaining the polarization degree and angle solely as functions of frequency.

The update of macroparticle spectra at shocks is based on the work of Vaidya et al. (2018) and Mukherjee et al. (2021). First, we implemented the shock-detection algorithm proposed in Mignone et al. (2012). We flag as shock the simulation cells where these two conditions are satisfied:

$$\nabla ·v<0\text{and}\frac{|\nabla P|}{P}>\zeta ,$$(35)

where ζ is an arbitrary threshold. When a macroparticle enters a flagged zone, the spectrum is not updated, and the pressure interpolated at each step is recorded. Once the macroparticle exits the flagged zone, the step with the lowest pressure is identified as “upstream”, while the step with the highest pressure is designated as “downstream”. The spectrum of the exiting macroparticle is then updated using the fluid quantities interpolated at upstream and downstream points. In Vaidya et al. (2018), the post-shock spectrum is reset to a power-law distribution, while its normalization and minimum energy are computed by requiring the macroparticle to have a fraction of the fluid energy and number density.

This approach presents two drawbacks: first, the macroparticle history, namely its spectrum temporal evolution before the shock, is forgotten, second, a computational cell could contain more than one shocked macroparticle and even non-shocked macroparticles; therefore, fluid energy and number density must be distributed among the shocked macroparticles taking into account these effects. We adopt the solution proposed in Mukherjee et al. (2021). Consider a particle entering the shock with spectrum ${\mathcal{N}}_{\mathrm{pre}}^{\prime }$, the post-shock spectrum ${\mathcal{N}}_{\mathrm{post}}^{\prime }$ is updated through a convolution

$${\mathcal{N}}_{\text{post}}^{\prime }(E^{\prime })=C{\int }_{E_{\text{min}}^{\prime }}^{E^{\prime }}{\mathcal{N}}_{\text{pre}}^{\prime }(E^{″}){\left(\frac{E^{\prime }}{E^{″}}\right)}^{-q+2}\frac{\mathrm{d}{E}^{″}}{E^{″}}·$$(36)

The upstream and downstream conditions determine the index q (see later), while the minimum energy $E_{\min }^{\prime }$ is fixed to the minimum energy of the pre-shock spectrum. With the convolution, the macroparticle history is not neglected. Specifically, in case a macroparticle encounters a strong shock followed by a weak shock, the spectrum is not anomalously steepened.

The normalization constant C is determined following a two-step procedure. First, we calculate the internal energy density to distribute among the shocked macroparticles

$$\mathrm{\Delta }E=f_E\rho ϵ-\sum _{\text{nsh}}\varepsilon ,$$(37)

where ρϵ is the fluid internal energy density, f_E is the fraction of the internal energy density available to the shocked macroparticles, and ∑_nshɛ is the sum of the energy densities of the non-shocked macroparticles. Note that ϵ is the fluid internal energy per unit of mass and it is determined by the Taub-Matthews equation of state since we are considering a relativistic fluid (Mignone et al. 2007). If ΔE<0, the shocked macroparticle spectra are not updated. Second, we compute the fluid number density to distribute among the shocked macroparticles

$$\mathrm{\Delta }N=f_N\frac{\rho }{\mu m_a}-\sum _{\mathrm{sh}\hspace{0.5em}\&\hspace{0.5em}\mathrm{nsh}}{N}_i,$$(38)

where m_a is the atomic mass unit, μ is the mean molecular weight for ionized gas, f_N is the fraction of the fluid number density available to the shocked macroparticles, and ∑_{sh & nsh}N_i is the sum of the number densities of non-shocked macroparticles along with the contribution of the shocked macroparticles, normalized in the previous step. If ΔN<0, the normalization constants C are reduced to obtain ΔN = 0.

3.1. Upstream and wall frame

In our framework, we aimed to refine the shock acceleration approach proposed in Vaidya et al. (2018). We calculated the shock quantities at MHD scales (e.g., the upstream and downstream magnetic fields) to derive quantities (e.g., upstream magnetization) that can be directly compared with PIC simulation results. In their turn, kinetic simulations provide key features of the post-shock spectrum, including the spectral index q, the maximum Lorentz factor γ_max, and the fractions of internal energy and particle number density available to shocked macroparticles, denoted as f_E and f_N, respectively.

However, not all shock configurations lead to efficient nonthermal particle acceleration. Given moderate magnetization and a relativistic shock, the upstream magnetic obliquity θ_up (i.e., the angle between the upstream magnetic field and the shock normal) plays a key role. In fact, in a magnetized medium, gyrating particles are constrained to move along the magnetic field lines, which in turn are advected downstream by the flow. Therefore, if the upstream magnetic obliquity is too large, the particles should move faster than the speed of light to return upstream (Begelman & Kirk 1990). This scenario is called superluminal shock, and kinetic simulations confirm that in this case, no acceleration takes place (Crumley et al. 2019 for the transrelativistic case, Sironi & Spitkovsky 2009, 2011 for the relativistic case). More quantitatively, a shock is defined as superluminal when cos θ_up<β_sh,up, where β_sh,up is the shock speed normalized to c (Kirk & Heavens 1989; Begelman & Kirk 1990). Both quantities are defined in the frame where the upstream is at rest. Only subluminal shocks can efficiently accelerate particles. PIC simulations show that in this scenario the post-shock spectrum displays a nonthermal power-law tail (Crumley et al. 2019 for the transrelativistic case, Sironi & Spitkovsky 2009, 2011 for the relativistic case). Ultimately, when a macroparticle crosses a shock, it is essential to compute the upstream magnetic inclination and shock speed. We need to transform the fluid quantities from the laboratory frame to the upstream rest frame, where the superluminal condition is defined.

Moreover, PIC simulations show that the features of the post-shock spectrum depend on the upstream magnetization, magnetic inclination, and Lorentz factor (Crumley et al. 2019 for the transrelativistic case, Sironi & Spitkovsky 2009, 2011 for the relativistic case). However, in PIC simulations the fluid quantities are computed in the wall frame, defined as the frame in which the normal component of the downstream velocity is null. Therefore, it is also necessary to transform the fluid quantities from the laboratory to the wall frame to take advantage of the results of the PIC simulations. For the implementation of the transformations from the laboratory to the upstream and wall frames, refer to Appendix C.

3.2. Cooling versus advection time

Before showing our results, we first checked the consistency of our framework. We are implicitly assuming that the shock structures observed in PIC simulations are much smaller than the scales we consider. For fluid quantities, this assumption is naturally satisfied, as their temporal and spatial evolution is described by the RMHD equations, which, by definition, treat the plasma on scales much larger than shock kinetic scales. Conversely, the situation is less clear for the scales governing the evolution of the macroparticle spectrum. In particular, the synchrotron cooling length of the nonthermal particle can be comparable with the decay length of the downstream microturbulence. In this scenario, acceleration and cooling would not be fully spatially separable since the contribution of the downstream region close to the shock surface to the cooling process is not negligible. A more quantitative evaluation is thus necessary (see Vanthieghem et al. 2020 for a similar discussion).

PIC simulations show that the downstream microturbulence decay length ℓ is on the order of 10³c/ω_p,i, where ω_p,i is the ion plasma skin depth. From fluid simulations, the thermal proton number density inside the jet is equal to n_p,th = 1 cm⁻³. Hence, we obtain ℓ∼10¹⁰cm. Assuming a strong shock, the advection velocity in the downstream region is v_adv≈c/3, so the time a particle spends in the microturbulence before being advected in the far downstream is t_adv∼3ℓ/c. The synchrotron cooling time t_cool immediately after the shock depends on the microturbulence magnetic field strength δB. From our fluid simulations, we can estimate that the downstream magnetic field in the fluid frame is $B_2^{\prime }\sim 5\times {10}^{-2}\mathrm{G}$. However, PIC simulations show that the ratio between the microturbulence and far downstream magnetic energy density can be significantly large (∼50, e.g., Tavecchio et al. 2018), which implies that δB∼10⁻¹ G. The comparison between the two timescales gives an upper limit for the Lorentz factor below which the microturbulence contribution to cooling is negligible:

$$\gamma <\frac{2\pi m_e{c}^2}{{\sigma }_T{B}^2\ell }\approx 7.7\times {10}^{10}{\left(\frac{\delta B}{{10}^{-1}\mathrm{G}}\right)}^{-2}{\left(\frac{\ell }{{10}^{10}\mathrm{cm}}\right)}^{-1}.$$(39)

The hypothesis that the cooling in the microturbulence region is negligible is valid in the Lorentz factor range of interest. Assuming the magnetic field in the emitting region $B_2^{\prime }\sim 5\times {10}^{-2}\mathrm{G}$ and a relativistic Doppler factor D∼10, the Lorentz factor of the electrons emitting in the keV emission or lower is ≲4×10⁵, well below the limit given by Eq. (39).

4. Results

In this section, we report the results of our post-processing framework, applied to simulations of recollimated magnetized jets.

Our post-processing procedure was applied to the axisymmetric stationary state, obtained from 2D simulations and interpolated to 3D. The steady state corresponds to a jet collimated by a higher-pressured environment, as shown in Fig. 1. Since the collimated jet is supersonic, a recollimation shock forms (Komissarov & Falle 1997). After this shock, the pressure increases, reaching an equilibrium with the external medium. Meanwhile, the Lorentz factor decreases, but because the recollimation shock is oblique, the flow remains relativistic. As streamlines converge toward the jet axis, they are reflected by the reflection shock. Because the incoming flow is nearly perpendicular to the reflection shock, the dissipation of kinetic energy is more significant at this stage. This results in a further pressure increase, causing the downstream region of the reflection shock to become over-pressured compared to the post-recollimation shock region and the external medium. Consequently, the jet re-expands. However, if the external medium maintains a stationary pressure profile, the jet undergoes another recollimation, leading to an oscillatory behavior of the jet radius. In this process, the relativistic fluid repeatedly experiences cycles of recollimation and reflection shocks. Bulk kinetic energy is cyclically converted to thermal energy through shock heating and then reconverted to kinetic energy via adiabatic cooling (Matsumoto et al. 2012). The confinement of jets, leading to the development of the chain of recollimation-reflection shocks, has been extensively studied using 2D relativistic hydrodynamic (RHD) simulations (e.g., Gomez et al. 1995; Agudo et al. 2001; Bodo & Tavecchio 2018) and as well as a few 2D RMHD simulations (Mizuno et al. 2015; Martí et al. 2016), which (at low-magnetization, as in our case) provide similar results to RHD ones. Our 2D RMHD simulations are in agreement with previous results. However, jets often exhibit more complex behavior due to various instabilities that depend on the jet intrinsic properties (Birkinshaw 1996; Bodo et al. 2013, 2019; Gourgouliatos & Komissarov 2018) and its surrounding environment (Porth & Komissarov 2015; Tchekhovskoy & Bromberg 2016). A full 3D simulation is required to study recollimation instabilities in detail, but such an analysis is beyond the scope of this paper (Boula et al., in prep.).

Fig. 1. Map of the pressure (left) and Lorentz factor (right) in the xz plane (Simulation A). The main structures are indicated. The cells tagged as shock regions in the xz plane are highlighted in light blue.

Below, we outline the selected values for the free parameters used in the post-processing code. We fixed ζ = 1 (see Eq. (35)), which is smaller than previously used values (for example, Mignone et al. 2007). After some tests, we observed that large values of ζ produce thin shock regions, especially at the reflection shock. By selecting a smaller value, we obtained thicker shock regions that allow a more accurate determination of upstream and downstream states, leading to improved estimates of the derived quantities, such as the shock normal and speed.

We must specify the line of sight direction, which is crucial for computing the relativistic Doppler factor. The viewing direction is defined by the polar θ (also called in the text observer viewing angle) and azimuthal angle ϕ. Due to the axisymmetric nature of our solutions, the azimuthal angle is irrelevant, so, without loss of generality, we set ϕ = 0. As for the polar angle, we fixed θ = 5°, suitable for blazars. For a discussion on the dependence of the polarimetric properties on the observing viewing angle, see Section 4.2.

Macroparticles were injected into the jet base at the start of the post-processing. Since we were considering a steady-state scenario, continuously injecting macroparticles over multiple simulation time steps is equivalent to an instantaneous injection at a single simulation time step. In this context, the full history of a macroparticle moving along the jet provides for macroparticles injected at later times. Two macroparticles were injected into each jet base cell and we selected the time step to ensure two macroparticles for every step in the z-direction, namely Δt = 0.005 z₀/c in code units. We fixed the maximum time t_max = 550Δt = 2.75 z₀/c. This ensures that the motion of most of the macroparticles ends before reaching the second recollimation shock, which was intentionally excluded from our analysis (see later).

Each macroparticle spectrum was sampled on an energy grid spanning from E_min = 3m_ec² to E_max = 10⁹m_ec². The energy grid is logarithmically spaced, with a resolution of 5 points per decade. Regarding the initial energy distribution of macroparticles, we opted for a steep power law $\mathcal{N}\propto E^{-\alpha }$ with α = 9, and limited the Lorentz factor range to γ∈[10¹,10²]. With this injection spectrum, the macroparticle history before the first shock is irrelevant. Even if shock conditions are not optimal for strong acceleration, the convolution in Eq. (36) assures that a slightly flatter index dominates, resulting in a power law with the just-mentioned index, effectively erasing the influence of the steep injection spectrum. The values chosen for the free parameters are summarized in Table 1.

We did not consider the entire simulation domain for our post-processing analysis. Instead, we focused on a reduced region defined by [−0.2,0.2]×[−0.2,0.2]×[1,4.5], in units of z₀. This choice allowed us to exclude most of the external medium, where no shocks, and therefore acceleration, occur. Moreover, we were exclusively interested in the first couple of recollimation-reflection shocks. As shown by Costa et al. (2024) in unmagnetized recollimated jets, the turbulence generated downstream of the reflection shock dissipates the bulk kinetic energy of the flow and it weakens the subsequent shocks, to finally suppress their formation. Moreover, we anticipated that the reflection shock is the site where most high-energy radiation is produced (see later). Since we expected a similar behavior for jets with low magnetization, we left out shocks beyond the first reflection shock in our analysis.

In Fig. 2, we present a 2D slice at y = 0 showing the Lorentz factor distribution of Simulation C. Overlaid on this, we display the trajectories of three randomly selected macroparticles. These macroparticles originate at the jet base and, after a few steps, encounter the first recollimation shock. Since the macroparticles move in three dimensions, their exact time steps of being tagged as “shocked” do not necessarily align with the shock region in the xz plane. However, if we imagine rotating the xz-shock structure around the z-axis, we obtain the full 3D shock region. This visualization clarifies that the macroparticles do indeed traverse the shock. Beyond the recollimation shock, the jet narrows, causing the macroparticles’ trajectories to shift closer to the z-axis as they follow the jet streamlines. Upon reaching the first reflection shock, their paths deviate outward again due to the jet's re-expansion. Finally, only the central macroparticle interacts with the second recollimation shock. As already pointed out, the post-processing duration was selected to ensure that only a small fraction of macroparticles reach the second recollimation shock.

Fig. 2. Map of the Lorentz factor in the xz plane with the projection of the trajectory of 3 macroparticles (number 799, 3983, 7789, Simulation C). The trajectory color represents the macroparticle status: white when inside a shock region and black when outside. The cells tagged as shock regions in the xz plane are highlighted in light blue.

4.1. Subluminal versus superluminal scenario

Every time a macroparticle crossed a shock, the first step was to identify whether the configuration is subliminal or superluminal. In the superluminal scenario, no particle acceleration occurs, and therefore the post-shock spectrum remains unchanged. Moreover, we computed the wall frame parameters required for spectrum updates, specifically the upstream magnetic field inclination, magnetization, and Lorentz factor, so that these values are readily available if needed. By comparing these quantities with PIC simulations, we could determine the slope of the post-shock energy distribution for subluminal cases. Figure 3 illustrates the distribution of the essential parameters used to classify shock crossings in Simulation A (refer to Appendix D for Simulation B and C). In addition, we plot each shock crossing in the parameter space defined by cos θ_up and β_sh,up. This allows us to visually distinguish which configurations are superluminal or subluminal.

Fig. 3. Left and middle panels: Distribution of the upstream magnetic inclination and shock speed in the upstream rest frame for Simulation A. Right panel: Comparison of these two quantities for each configuration. The dot color indicates the z coordinate of the macroparticle exiting point from the shock. The configurations inside the red region are superluminal.

Except for a few outliers, likely due to inaccuracies in reconstructing the upstream rest frame parameters, almost all configurations are superluminal in all scenarios. For Simulation A and B, at low z (blue points) recollimation shock crossings lie deep within the superluminal region, and as the crossing occurs at higher z (green points), the parameters gradually move closer to the boundary between superluminal and subluminal regimes. Near the top of the recollimation shock, some configurations even become subluminal in Simulation A (dark green points). This behavior can be explained by the relative inclination between the shock normal and the streamlines. The inclination is high for low z, resulting in full superluminal configurations. As z grows, the inclination decreases, leading to shock crossing near the superluminal-subluminal limit. The reflection shock behaves oppositely. At low z, the shock crossings remain close to the superluminal-subluminal boundary, with a few subluminal configurations in Simulation A (dark red points), but as z increases, the shock crossings revert to fully superluminal conditions (red points). Overall, the recollimation shock transitions from strongly superluminal behavior toward the subluminal threshold with increasing z, while the reflection shock moves from near-sub to firmly superluminal conditions as z grows. Compared to the recollimation shock, the reflection shock speed is lower, a trend also evident for the few macroparticles that encounter the second recollimation shock (yellow points). In Simulation C, all configurations exhibit a strong superluminal behavior. This is primarily due to the magnetic field configuration, as the field is predominantly toroidal. Consequently, the upstream magnetic inclination remains low in nearly all cases.

Therefore, we conclude that, assuming a laminar flow, no particle acceleration can occur in the shocks associated with recollimation.

4.2. Turning on acceleration

No acceleration is possible, except for very few macroparticles. Thus, we decided to change our starting assumptions. Previously, we were supposing that the fluid is laminar on all lengths, from fluid to kinetic scales. The distinction between superluminal and subluminal configuration is crucial if the upstream flow is laminar (Crumley et al. 2019 for the transrelativistic case, Sironi & Spitkovsky 2009, 2011 for the relativistic case). Recent simulations indicate that the presence of turbulence or inhomogeneities on kinetic scales in the upstream medium can reignite acceleration, even in magnetized and superluminal configurations. In one scenario, turbulence frees particles from tightly following magnetic field lines, preventing their advection into the downstream region and thus enabling acceleration at the shock (Bresci et al. 2023). In another scenario, the shock interaction with these inhomogeneities favors the formation of intense downstream turbulence, which accelerates particles, leading to the formation of a power-law tail in the post-shock spectra (Demidem et al. 2023). For both scenarios, the parameter space is not extensively explored, as for laminar flows. Therefore, we could not compare the shock crossing parameters directly with the simulations to deduce the post-shock properties. For this reason, we were forced to use phenomenological formulae.

For the power law index, we opted for a generalization of the classical formula (Drury 1983), which was originally derived under the assumption of isotropic diffusion in the particle velocity angle relative to the shock normal (Keshet & Waxman 2005),

$$q=\frac{3{\beta }_{1,\mathrm{NIF}}-2{\beta }_{1,\mathrm{NIF}}{\beta }_{2,\mathrm{NIF}}^2+{\beta }_{2,\mathrm{NIF}}^3}{{\beta }_{1,\mathrm{NIF}}-{\beta }_{2,\mathrm{NIF}}}-2,$$(40)

where β_1,NIF and β_2,NIF are respectively the upstream and downstream velocity normalized to c in the normal incident frame (NIF), where the shock is at rest and the upstream velocity is directed along the shock normal (for more details about the transformation from the laboratory frame to the NIF, see Appendix C). This formula is justified if strong upstream turbulence detaches particles from magnetic field lines, like in the scenario investigated in Bresci et al. (2023), effectively making the shock behave as if it were low magnetized, as assumed in Keshet & Waxman (2005).

Regarding the maximum energy of the post-shock spectra, we followed the standard approach (e.g., Kirk et al. 1994; Böttcher & Dermer 2010; Baring et al. 2017). Supposing that the acceleration time is a multiple of the gyrofrequency (i.e., t_a=λ2πr_g/c; this corresponds to assume a “gyro-Bohm” spatial diffusion coefficient, e.g., Kirk et al. 1994) and the main cooling process is synchrotron emission, the maximum Lorentz factor, dictated by equilibrium between energy gains and losses, is given by

$${\gamma }_{\max }=\frac{E_{\max }}{m_e{c}^2}={\left(\frac{9m_e^2{c}^4}{8\pi B^{\prime }\lambda e^3}\right)}^{\frac{1}{2}},$$(41)

where B′ is the fluid rest frame magnetic field interpolated at the point where the macroparticle exits the shock region. The free parameter λ, also called acceleration efficiency, contains the shock microphysics and can be interpreted as the ratio between the diffusion mean free path and the particle gyroradius. We set λ = 1000, while in Vaidya et al. (2018) this parameter is derived through a more complex approach, which gives, on average, λ∼1, namely the so-called Bohm limit. As shown in Garson et al. (2010), Boettcher et al. (2012), Summerlin & Baring (2012), this value implies extremely efficient shocks, resulting in electrons producing synchrotron radiation up to the MeV band. Thus, we adopted a significantly larger λ, ensuring that synchrotron emission peaks inside the X-ray band, as observed for HBLs and EHBLs (in Appendix E we report the flux as a function of frequency for the three cases). On kinetic scales, this value corresponds to relatively weak microturbulence around the shock, resulting in a slower diffusion and, thus, a higher acceleration time. This is in tension with the previously described scenario, where particles are either detached by magnetic field lines or directly accelerated by strong turbulence. However, the maximum particle energy could still be limited by the spatial scale of the microturbulence. As shown in Sironi et al. (2013), if a particle gyroradius exceeds the scale of the microturbulence, it will be advected away from the shock front, effectively halting the acceleration process¹. This mechanism could alleviate the discrepancy between the required acceleration efficiency for effective particle acceleration and the observed synchrotron emission limits. While this effect has been verified in PIC simulations of laminar flows, it has not yet been fully explored in turbulent or inhomogeneous flow conditions; therefore, we stick to Eq. (41).

For a laminar flow, the ratio of nonthermal to thermal energy density is 10⁻³ for a transrelativistic shock (Crumley et al. 2019) and 10⁻¹ for a relativistic shock (Sironi & Spitkovsky 2009, 2011). The ratio of nonthermal to thermal number density is only available for the relativistic scenario, where it is 5×10⁻² (Sironi & Spitkovsky 2009, 2011). In the case of turbulent and inhomogeneous scenarios, these values are unavailable due to the limited runtime of the simulations. Given this, we adopt intermediate values of f_N=f_E = 10⁻². However, it is important to note that since we are primarily concerned with intensity ratios, these specific values do not affect the final results.

Figs. 4 and 5 represent a macroparticle trajectory and its spectral time evolution. The macroparticle is injected into the jet base, then it is advected by the flow until it reaches the recollimation shock. Before that, the spectrum is not visible in Fig. 5 because of the low normalization. Thanks to shock acceleration, the macroparticle spectrum becomes a power law, but, as it is advected from the recollimation to the reflection shock, an exponential cut-off emerges for high Lorentz factors due to the synchrotron cooling. When the macroparticle encounters the reflection shock, it is accelerated again. Since the magnetic field is higher due to the double shock compression, the maximum Lorentz factor is lower. On the other side, the internal energy and particle number density increase after both shocks, resulting in a larger normalization for the post-reflection shock spectrum. Finally, the macroparticle is again advected by the flow and cooled down due to synchrotron and, to a smaller extent, adiabatic losses.

Fig. 4. Map of the Lorentz factor in the xz plane with the projection of the trajectory of a macroparticle (number 1135, Simulation B). The color of the trajectory represents the time step, see Fig. 5.

Fig. 5. Time evolution of the normalized spectrum of the macroparticle in Fig. 4.

Figs. 6, 7, and 8 show the emissivity in the observer frame for Simulation B of a slice centered on the xz plane in the radio, optical, and X-ray band. The emissivities of Simulation A and C present similar features (see Appendices F and G). The evident asymmetry of the emissivity is due to the combined effect of relativistic beaming and the $|B^{\prime }\times {\hat{\text{n}}}^{\prime }|$ term in the synchrotron formulae, see Eqs. (18) and (19). The emission from portions of the flow with streamlines closely oriented toward the observer (e.g., the left side of the recollimation shock downstream) is more amplified, compared to the emission associated with the plasma whose velocity is oriented in other directions.

We neglected Faraday rotation and synchrotron self-absorption, which are significant below 10¹¹ Hz (Tavecchio et al. 2010b). Consequently, our model is not suitable for accurately reproducing radio emission from regions at parsec-scale distances. However, the observed radio emission is primarily produced by regions located farther from the supermassive black hole.

Fig. 9 displays the polarization parameters of the three configurations. In all cases, we obtain a strong chromaticity. In Simulations A and B, the polarization degree increases with the frequency, as observed by multiwavelength polarization campaigns of HBLs and EHBLs in a quiescent state. In particular, the polarization degree values of Simulation B are comparable with the IXPE and multiwavelength results. The chromaticity can be traced back to the structure of the synchrotron emissivity and the magnetic fields. In the radio band (Fig. 6), the region where the emissivity is significant (from now on emission region) includes both the downstream regions of the recollimation and the reflection shocks. In these regions, both components of the field, toroidal and poloidal, are significant (Fig. 10), thus reducing the total polarized fraction of the emission. On the other hand, the contribution of the recollimation shock downstream is less important in the optical band (Fig. 7), which in turn increases the degree of polarization since the macroparticles experience a more coherent magnetic field. For the same reason, the X-ray polarization degree is even higher, since the emission region shrinks further (Fig. 8), resulting almost exclusively in the downstream of the reflection shock. We recover here the main ingredient of the stratified shock scenario.

Fig. 6. Normalized observer frame synchrotron emissivity in the radio band of a 10-code units slab centered on xz plane (Simulation B).

Fig. 7. Same as Fig. 6 but in the optical band.

Fig. 8. Same as Fig. 6 but in the X-ray band.

The emissivity of the reflection shock downstream prevails for higher frequencies for two reasons. First, since the upstream converging flow is almost always perpendicular to the reflection shock, the kinetic energy dissipation is high, with a large amount of internal energy available for the accelerated macroparticles. Second, the magnetic field downstream of the recollimation shock is larger than in the recollimation post-shock region, thus the macroparticles have a larger synchrotron output.

Simulation C displays the opposite trend, namely the polarization degree decreases with frequency. In this case, the downstream of both the recollimation and reflection shock is dominated by the toroidal component; therefore, the polarization depends exclusively on the asymmetry of the emission region. In fact, since we are observing at a non-negligible angle, relativistic beaming and the $|B^{\prime }\times {\hat{\text{n}}}^{\prime }|$ term break the jet symmetry. As in the other two simulations, the downstream regions of both the recollimation and reflection shocks contribute to the radio radiative output. Since the radio emission region extends far from the jet axis, the asymmetry becomes more pronounced, leading to a higher degree of polarization. In contrast, the emission region in the optical and X-ray bands is more confined, eventually shrinking to just the downstream region of the reflection shock. Because asymmetry is less significant closer to the jet axis, the overall polarization degree is reduced.

In Fig. 9, the polarization angle is measured clockwise from the projection of the jet axis onto the sky plane. In Simulation C, the toroidal component of the fluid rest frame magnetic field is dominant by default. As a result, the polarization angle is aligned with the sky projection of the jet axis across all frequencies. In contrast, in Simulation A, the poloidal component of the magnetic field dominates. This causes the polarization angle to be perpendicular to the projection of the jet axis onto the sky plane.

Fig. 9. Polarization degree (top panel) and angle (bottom panel) of Simulation A (blue circles), B (orange triangles), and C (green squares). The polarization angles are measured clockwise from the projection of the jet axis onto the sky plane (solid gray line).

Fig. 10. Map of the poloidal (left panel) and toroidal (right panel) component of the fluid rest frame magnetic field in the xz plane (Simulation B).

In Simulation B, at high frequencies, the polarization angle is perpendicular to the jet axis projection, while at lower frequencies it is aligned with it. During the transition, the polarization degree temporarily vanishes before increasing again at even lower frequencies. This smooth transition in Simulation B contrasts with the abrupt changes observed in simpler models, such as cylindrical jets (e.g., Lyutikov et al. 2005). Due to the complex structure of the jet flow and the magnetic field in our simulation, the polarization angle can assume intermediate values. The polarization properties of Simulation B can be explained by the structure of the synchrotron emissivity and magnetic fields. The emission region at high frequencies corresponds primarily to the reflection shock downstream, where the poloidal component is significant. At lower frequencies, however, the emission is increasingly influenced by the recollimation shock downstream, where the toroidal component prevails. As a result, the growing contribution of the toroidal component reduces the polarization degree and causes the polarization angle to rotate. Finally, since the toroidal contribution is dominant for even lower frequencies, the polarization angle aligns with the jet axis projection and the polarization degree returns to increase.

Finally, in Appendix H, we report for Simulation A the dependence of the polarization degree as a function of the frequency on the observer viewing angle. The polarization degree varies smoothly across different frequencies as the viewing angle changes. At a fixed frequency, the polarization degree decreases as the viewing angle becomes smaller. At lower viewing angles, the emission region appears more symmetric because the beaming is more uniform. Since the magnetic field is axisymmetric, this results in a reduced polarization degree. However, as the viewing angle changes, the chromaticity between remains almost constant.

5. Discussion

In this paper, we used the Lagrangian macroparticle approach to probe the nonthermal emission of weakly magnetized blazar jets, specifically to replicate and interpret multifrequency polarimetric observations of HBLs and EHBLs. The Lagrangian approach bridges the gap between fluid-scale structures and kinetic-scale phenomena. While fluid simulations provide bulk velocities, fields, and shock structures at large scales, the macroparticle approach incorporates shock acceleration.

The first result we found is that shocks from a recollimating jet are characterized by superluminal configurations, limiting the effectiveness of acceleration in purely laminar flows. However, the study of alternative scenarios (explored through PIC simulations) that include turbulence and/or inhomogeneities on kinetic scales demonstrates that acceleration can occur, resulting in nonthermal power law distributions (e.g., Bresci et al. 2023; Demidem et al. 2023).

Assuming that acceleration is active, our results for Case B, for which the poloidal and the toroidal components are comparable, prove that the model can easily reproduce the chromaticity of the polarization degree observed in HBLs and EHBLs in quiescent states (e.g., Liodakis et al. 2022; Ehlert et al. 2023), even assuming globally ordered fields (i.e., no turbulence at fluid scales). This is in line with the findings of Bolis et al. (2024). Moreover, the model also naturally accounts for the observed overlapping of the optical and X-ray EVPA.

The model does not reproduce the alignment of the polarization angle with the jet axis observed in some HBL (e.g., Liodakis et al. 2022). This misalignment can be traced to the intrinsic magnetic field structure in our simulations. The optical and X-ray emission regions almost coincide with the reflection shock downstream. Specifically, the emission is higher around the z axis, where the toroidal magnetic field component is zero by definition, leaving the poloidal component dominant at these frequencies. As a result, the polarization angle is oriented perpendicular to the projected jet axis. A jet with a dominant toroidal component (case C) displays the EVPA aligned with the projected jet axis, but the properties of the degree of polarization (decreasing with increasing frequencies) are incompatible with the observations. It is worth noting here that the alignment is evident for Mrk501 (Liodakis et al. 2022), while it appears less definitive for Mrk421 (Di Gesu et al. 2022) and 1ES 0229+200 (Ehlert et al. 2023). At any rate, the mismatch could be explained by the jet bending between the inner X-ray and optical emission region and the outer extended region observed by radio telescopes (Di Gesu et al. 2022). Moreover, jets can change direction in time (Kostrichkin et al. 2025), thus requiring simultaneous radio imaging for an appropriate estimate of the alignment of the polarization angles with the jet projection.

Our results could depend on the magnetic field geometry since different profiles (e.g., Hu et al. 2025) could lead to subluminal shocks or different polarization features. We plan to explore different configurations in a future paper.

Our current post-processing code has so far been applied only to axisymmetric RMHD simulations without accounting for the system time evolution after performing 3D interpolation. As previously discussed, fluid-level turbulence can arise downstream of the reflection shock. While particle acceleration at shocks should not be affected, since in fully time-evolved 3D simulations the reflection shock remains largely intact (Costa et al. 2024; Boula et al., in prep.), the turbulent magnetic field of the reflection shock downstream could affect the emission. This could help to relax the fine-tuning of the pitch angle needed to reproduce observations. While Simulation B closely reproduces the IXPE and multiwavelength data, in Simulation A the polarization degree is too large and not enough chromatic. However, a more disordered magnetic field would reduce the overall polarization degree. Moreover, turbulence should have a greater impact on electrons that travel longer distances (namely emitting at lower frequencies), accentuating the chromaticity. Finally, although our present code only explains the polarization of blazars in a quiescent state, incorporating the time evolution and subsequent turbulence-driven symmetry breaking could help explain additional phenomena, such as the observed X-ray polarization angle swing in Mrk 421. Finally, turbulence may also contribute to particle acceleration. The current spectral update implementation is already designed to incorporate stochastic acceleration by turbulence (see also Kundu et al. 2021).

Another source of turbulence could be intrinsic to the jet flow. At fluid scales, our RMHD simulations display a laminar flow. However, the pre-shock fluid could be inhomogeneous and/or turbulent even on fluid scales. The non-laminarity of the fluid could result in a more turbulent downstream, leading also to magnetic field amplification (Giacalone & Jokipii 2007; Mizuno et al. 2014). As a result, the emission from the macroparticles could change, especially the predicted polarization degree. Moreover, strong turbulence could also further accelerate particles, as outlined in the previous paragraph.

It is also important to note that the framework we have described is well-suited for HBLs and EHBLs. For other classes of blazars, such as FSRQs, LBLs, and IBLs, the situation is more complex. These sources tend to have higher luminosities and lower characteristic peak frequencies, which implies that the conditions in the emitting regions differ significantly. They require larger magnetic fields and a higher nonthermal electron number density. Under such conditions, the maximum Lorentz factor for which the cooling length is higher than the decay length of the microturbulence, see Eq. (39), could be much lower, making our approximation inconsistent in the bands of interest. Moreover, the ratio of the high and low energy peaks is higher than for HBLs and EHBLs, thus the contribution of the inverse Compton cooling cannot be neglected for these sources. We plan to implement inverse Compton cooling and emission in future work (for a related example of external Compton processes on the cosmic microwave background, see Vaidya et al. 2018).

Acknowledgments

We thank E. Sobacchi for fruitful discussions. We acknowledge financial support from an INAF Theory Grant 2022 (PI F. Tavecchio). This work has been funded by the European Union-Next Generation EU, PRIN 2022 RFF M4C21.1 (2022C9TNNX). This work has been funded by ASI under contract 2024-11-HH.0. We acknowledge support by CINECA, through ISCRA and Accordo Quadro INAF-CINECA, and by PLEIADI, INAF – USC VIII, for the availability of HPC resources. Plots have been made with Python. We have used the following libraries: Numpy (Harris et al. 2020), Matplotlib (Hunter 2007), Scipy (Virtanen et al. 2020), and PyPluto (Mattia et al. 2025).

References

Appendix A: Derivation of the magnetic field profile

In this work, we simulated a helical magnetic field in a conical jet. In such a case, the field is made both of a toroidal field, perpendicular to the streamlines, and of a poloidal component, directed along the spherical radius r; no component is null in cylindrical coordinates: $\text{B}=\left(B_R,B_z,B_{\varphi }\right)$. We remind that a factor $\sqrt{4\pi }$ is reabsorbed in the definition of the magnetic field and that we work with c = 1. To find a simple way to prescribe analytically the magnetic field, it is easier to work in spherical coordinates (r,θ,ϕ), because $\text{v}=\left(v_r,0,0\right).$ We can assume that B_θ = 0, so that the field is $\text{B}=\left(B_r,0,B_{\varphi }\right)$; the poloidal and toroidal components are simply B_pol=B_r and B_tor=B_ϕ, and due to axisymmetry none of them depend on ϕ. These quantities are defined in the laboratory-frame; the proper frame magnetic field is B′=(B_r,0,B_ϕ/Γ). While the proper-frame electric field is null in ideal magneto-hydrodynamics, the laboratory frame electric field is given by

$$\text{E}=-\text{v}\times \text{B}=v_r{B}_{\varphi }{\hat{e}}_{\theta }.$$(A.1)

In ideal MHD, the magnetic flux is conserved. Since the jet's initial (and boundary) geometry is conical, its section is Σ∝r² and this determines a condition on the poloidal field

$${\int }_{r^2}\text{B}·d\text{S}=\text{conserved}\to B_r=\frac{g(\theta )}{r^2};$$(A.2)

on the other hand, the conservation of current implies

$$I={\oint }_{\partial \mathrm{\Sigma }}\text{B}·d\text{l}-\frac{d}{\mathit{dt}}{\int }_{\mathrm{\Sigma }}\text{E}·d\text{S}\propto {\oint }_r\text{B}·d\text{l}-k\frac{d}{\mathit{dt}}{\int }_{r^2}\text{E}·d\text{S},$$(A.3)

$$\to B_{\varphi }=\frac{h(\theta )}{r},$$(A.4)

where k is a constant and the displacement current is zero because E_r = 0. The dependence on the spherical radius is different for the two magnetic field components, so it is not possible to define a magnetic field in an equilibrium configuration in the whole conical jet. Since the initial configuration is merely a starting point, and the only region where the initial condition is maintained throughout the simulation is the nozzle, we limit to setting the equilibrium across the jet (in the θ direction) at its injection in (r = 1,θ):

$$\left[{\left[(\nabla ·\text{E})\text{E}+(\nabla \times \text{B})\times \text{B}\right]}^{\theta }-\frac{\partial p}{\partial \theta }\right](r=1,\theta )=0.$$(A.5)

Making explicit the pressure dependency on the spherical radius, we have $p=\left(\tilde{p}(\theta )+{\tilde{p}}_0\right)r^{-2\gamma }$, where ${\tilde{p}}_0$ is a constant, and γ_ad = 5/3 is the adiabatic index. Then, using also eq. A.1, we can solve the previous equation,

$$\frac{\partial }{\partial \theta }\left(\frac{h^2(\theta )}{{\mathrm{\Gamma }}^2}\right)+2cot\theta \left(\frac{h^2(\theta )}{{\mathrm{\Gamma }}^2}\right)=-\frac{\partial g^2(\theta )}{\partial \theta }-\frac{\partial \tilde{p}(\theta )}{\partial \theta },$$(A.6)

$$\frac{\partial }{\partial \theta }\left(\frac{h^2(\theta )}{{\mathrm{\Gamma }}^2}\right)+2cot\theta \left(\frac{h^2(\theta )}{{\mathrm{\Gamma }}^2}\right)=-\frac{\partial l^2(\theta )}{\partial \theta },$$(A.7)

where we defined for simplicity $l^2(\theta )=g^2(\theta )+\tilde{p}(\theta )$.

If we adopt the following Gaussian profile $l^2(\theta )=e^{-2{\mathcal{X}}^2}$, where $\mathcal{X}=\left(cos\theta -1\right)/{\psi }_{\chi }$ and ψ_χ is an arbitrary constant, equation A.7 is solved by

$$h^2(\theta )={\mathrm{\Gamma }}^2\left[-{\text{e}}^{-2{\mathcal{X}}^2}-\frac{{\psi }_{\chi }}{2{sin}^2\theta }\left({\psi }_{\chi }-{\psi }_{\chi }{\text{e}}^{-2{\mathcal{X}}^2}+2\sqrt{\pi }\text{erf}\left(\sqrt{2}\mathcal{X}\right)\right)\right].$$(A.8)

The functions l(θ) and h(θ) yield appropriate profiles for the poloidal and toroidal field components if ψ_χ = 10⁻²θ_j, as shown in Fig. A.1 for θ_j = 0.1: max[l(θ)]=l(0) = 1, h(0) = 0, max[h(θ)]≃h(0.05)≃0.7 Γ, and both functions decay to 0 for θ>θ_j = 0.1.

Fig. A.1. Profiles of l(θ) and h(θ)/Γ.

We define the lab-frame toroidal field as

$$B_{\varphi }=\frac{B_0\mathrm{\Gamma }}{r}\sqrt{-{\text{e}}^{-2{\mathcal{X}}^2}-\frac{{\psi }_{\chi }}{2{sin}^2\theta }\left({\psi }_{\chi }-{\psi }_{\chi }{\text{e}}^{-2{\mathcal{X}}^2}+\sqrt{2\pi }\text{erf}\left(\sqrt{2}\mathcal{X}\right)\right)},$$(A.9)

and the poloidal field as

$$B_r=\alpha B_0\frac{e^{-{\mathcal{X}}^2}}{r^2},$$(A.10)

with the parameter α that determines the pitch. The pressure profile is chosen to satisfy equation A.7,

$$\tilde{p}(\theta )=(1-{\alpha }^2)B_0^2{e}^{-2{\mathcal{X}}^2}\to p(r,\theta )=(1-{\alpha }^2)\frac{B_0^2}{2}\frac{{\text{e}}^{-2{\mathcal{X}}^2}}{r^{2{\gamma }_{\text{ad}}}}+\frac{p_{j,0}}{r^{2{\gamma }_{\text{ad}}}},$$(A.11)

where the constant p_j,0 is chosen as floor jet pressure value at injection.

In the previous calculations, we found a magnetic field configuration that is in equilibrium in the θ direction at (r = 1,θ). Instead, simulations are performed in a cylindrical domain, and the jet injection boundary corresponds to the region (R<θ_j,z = 1), where R=r sin θ and z=r cos θ. For small values of θ, within the jet, where the profile is relevant, we can approximate the two equilibrium states, since R(r = 1,θ)∼θ and z(r = 1,θ)∼1, so we inject the cylindrical field

$$\text{B}(R,z)=\left(B_r(R,z)\frac{R}{\sqrt{R^2+z^2}},B_{\varphi }(R,z),B_r(R,z)\frac{z}{\sqrt{R^2+z^2}}\right).$$(A.12)

Appendix B: Numerical scheme for spectral update

In our code, the numerical scheme for the spectral update step is a second-order implicit-explicit Runge-Kutta scheme, specifically the strong stability preserving algorithm (2,2,2), see Pareschi & Russo (2005). The choice of this scheme permits the future introduction of a diffusion term accounting for turbulence acceleration, as in Kundu et al. (2021). Since the advection term is solely present, the scheme reduces to

$$\{\begin{array}{cc}\hfill & {\chi }_{i,*}={\chi }_{i,n}+\mathrm{\Delta }{t}_n{\mathcal{A}}_{i,n},\hfill \\ \hfill & {\chi }_{i,n+1}={\chi }_{i,n}+\frac{\mathrm{\Delta }{t}_n}{2}\left[{\mathcal{A}}_{i,n}+{\mathcal{A}}_{i,*}\right].\hfill \end{array}$$(B.1)

The subscript i indicates the i-th point in the logarithmically equally spaced energy grid, while $\mathcal{A}$ is the advection term, discretized as

$${\mathcal{A}}_{i,n}=-{\xi }_i^{\prime }\frac{{\mathcal{F}}_{i+\frac{1}{2},n}^{\text{adv}}-{\mathcal{F}}_{i-\frac{1}{2},n}^{\text{adv}}}{\mathrm{\Delta }\xi },$$(B.2)

where ξ_i=log(E_i/E_min)/log(E_max/E_min) and ${\xi }_i^{\prime }=[E_{i}log(E_{\max }/E_{\min }){]}^{-1}$ is the corresponding Jacobian. ℱ^adv indicates the advection flux, discretized as

$${\mathcal{F}}_{i+\frac{1}{2},n}^{\text{adv}}=\{\begin{array}{cc}H({\gamma }_{i+\frac{1}{2}}){\chi }_{i+\frac{1}{2},n}^L\hfill & \text{if}\hspace{0.5em}H({\gamma }_{i+\frac{1}{2}})>0,\hfill \\ H({\gamma }_{i+\frac{1}{2}}){\chi }_{i+\frac{1}{2},n}^R\hfill & \text{if}\hspace{0.5em}H({\gamma }_{i+\frac{1}{2}})<0,\hfill \end{array}$$(B.3)

where H is the sum of the two terms in the round brackets of Eq. (15). Finally, the left χ^L and right χ^R are defined as

$${\chi }_{i+\frac{1}{2}}^L={\chi }_i+\frac{\delta {\chi }_i}{2},$$(B.4)

$${\chi }_{i+\frac{1}{2}}^R={\chi }_{i+1}-\frac{\delta {\chi }_{i+1}}{2},$$(B.5)

$$\text{with}\hspace{0.5em}\delta {\chi }_i=\{\begin{array}{cc}\frac{2\mathrm{\Delta }{\chi }_{i+\frac{1}{2}}\mathrm{\Delta }{\chi }_{i-\frac{1}{2}}}{\mathrm{\Delta }{\chi }_{i+\frac{1}{2}}+\mathrm{\Delta }{\chi }_{i-\frac{1}{2}}}\hfill & \text{if}\hspace{0.5em}\mathrm{\Delta }{\chi }_{i+\frac{1}{2}}\mathrm{\Delta }{\chi }_{i-\frac{1}{2}}>0\hfill \\ 0,\hfill & \text{otherwise},\hfill \end{array}$$(B.6)

where $\mathrm{\Delta }{\chi }_{i\pm \frac{1}{2}}=\pm \left({\chi }_{i\pm 1}-{\chi }_i\right)$. Note that for this scheme two ghost cells are needed on both sides of the energy grid. To conserve particle number, a zero flux condition must be imposed at both borders. For our algorithm, this translates into the following boundary conditions

$${\chi }_i^{\text{adv}}=\{\begin{array}{cc}-{\chi }_{2i_b-i-1}^{\text{adv}},\hfill & \text{for}\hspace{0.5em}i<i_b,\hfill \\ -{\chi }_{2i_e-i+1}^{\text{adv}},\hfill & \text{for}\hspace{0.5em}i>i_e,\hfill \end{array}$$(B.7)

$$\text{with}\hspace{0.5em}{\mathcal{F}}_{i_b-\frac{1}{2}}^{\text{adv}}={\mathcal{F}}_{i_e+\frac{1}{2}}^{\text{adv}}=0,$$(B.8)

where i_b and i_e are the indices of the first and last points of the energy grid, respectively. This scheme is second-order accurate, but the Courant–Friedrichs–Lewy condition ($C=\mathcal{A}\mathrm{\Delta }t/\mathrm{\Delta }E<1$, with C the Courant number) must be satisfied since the advection term is treated explicitly.

Appendix C: From laboratory to upstream and wall frame: detailed implementation

To boost to the wall and upstream rest frames, we first compute the shock normal in the laboratory frame ${\hat{\text{n}}}_{\mathrm{sh}}$ by means of the coplanarity theorem (it can be proved directly from the RMHD jump conditions), which states that the upstream and downstream magnetic field along with the velocity jump lie in the same plane as the shock normal (Landau & Lifshitz 1960),

$${\hat{\text{n}}}_{\mathrm{sh}}=\{\begin{array}{cc}\hfill & \pm \frac{{\beta }_2-{\beta }_1}{|{\beta }_2-{\beta }_1|}\text{parallel shocks},\hfill \\ \hfill & \pm \frac{(B_1\times \mathrm{\Delta }\beta )\times \mathrm{\Delta }B}{|(B_1\times \mathrm{\Delta }\beta )\times \mathrm{\Delta }B|}\text{otherwise},\hfill \end{array}$$(C.1)

where the subscripts 1 and 2 indicate, respectively, the upstream and downstream values. Note that a different routine is implemented for parallel shocks since ΔB = 0. Then, by enforcing number flux conservation across the shock surface, we calculate the shock speed in the laboratory frame, whose direction is along the shock normal:

$${\beta }_{\mathrm{sh}}=\frac{{\rho }_2{\gamma }_2{\beta }_2·{\hat{\text{n}}}_{\mathrm{sh}}-{\rho }_1{\gamma }_1{\beta }_1·{\hat{\text{n}}}_{\mathrm{sh}}}{{\rho }_2{\gamma }_2-{\rho }_1{\gamma }_1}.$$(C.2)

Given the shock normal and speed, the fluid quantities can be transformed to the normal incident frame (NIF), where the shock is at rest and the upstream velocity is directed along the shock normal. This transformation requires a two-step procedure, in analogy with the process outlined in Kirk & Heavens (1989) and Summerlin & Baring (2012) for the de Hoffmann-Teller frame. First, we perform a Lorentz boost equal to the shock speed and along the shock normal. Then, the second boost is along the transverse component of the upstream velocity. The two-step procedure has an important advantage: the shock normal does not vary across the three different frames, namely ${\hat{\text{n}}}_{\mathrm{sh}}={\hat{\text{n}}}_{\mathrm{sh},\mathrm{int}}={\hat{\text{n}}}_{\mathrm{sh},\mathrm{NIF}}$, where the subscript “int” indicates the intermediate frame between the laboratory and NIF. The first boost is along the shock normal, which therefore remains unchanged, while the second boost, though perpendicular to the normal, also leaves the normal unchanged, since the shock is at rest. Defining Λ[u] as a generic Lorentz boost along the velocity u, a four-vector x in the laboratory frame can be transformed to the NIF as follows

$${\text{x}}_{\mathrm{NIF}}=\mathrm{\Lambda }[v_{1,\mathrm{int}}-(v_{1,\mathrm{int}}·{\hat{\text{n}}}_{\mathrm{sh}}){\hat{\text{n}}}_{\mathrm{sh}}]\mathrm{\Lambda }[v_{\mathrm{sh}}{\hat{\text{n}}}_{\mathrm{sh}}]\text{x},$$(C.3)

We can transform to the upstream rest frame with a further step. We perform a boost along the NIF upstream velocity, which is parallel to the shock normal:

$${\text{x}}_{\mathrm{up}}=\mathrm{\Lambda }[v_{1,\mathrm{NIF}}]{\text{x}}_{\mathrm{NIF}}.$$(C.4)

The boost is again along the shock normal, which, therefore, remains unchanged. Hence, we can substitute ${\hat{\text{n}}}_{\mathrm{sh},\mathrm{up}}$ with ${\hat{\text{n}}}_{\mathrm{sh}}$ in the superluminal condition.

Similarly, it is possible to transform from the NIF to the wall frame in one step, with a single boost along the normal component of the downstream velocity:

$${\text{x}}_{\mathrm{wall}}=\mathrm{\Lambda }[(v_{2,\mathrm{NIF}}·{\hat{\text{n}}}_{\mathrm{sh}}){\hat{\text{n}}}_{\mathrm{sh}}]{\text{x}}_{\mathrm{NIF}}.$$(C.5)

As before, the shock normal does not vary, namely ${\hat{\text{n}}}_{\mathrm{sh},\mathrm{wall}}={\hat{\text{n}}}_{\mathrm{sh}}$. Moreover, the upstream velocity remains along the shock normal. Note that, since the magnetization is low, it is possible to prove that the transverse component of the downstream velocity is nonrelativistic; therefore, no extra boost is needed to move from the wall frame to the downstream rest frame.

Transforming from the laboratory to the NIF, upstream, and downstream rest frame may seem unnecessarily complex. A single Lorentz boost appears to be the simplest approach for the upstream and downstream rest frames. However, this approach results in a rotation of the shock normal (Kirk & Heavens 1989; Ballard & Heavens 1991). While this difference in the shock normal behavior may seem like a discrepancy, it is actually a direct consequence of the well-known Thomas-Wigner rotation: the composition of two non-collinear boosts results in a Lorentz transformation that is not a pure boost but is the composition of a boost and a rotation. To avoid complications arising from this effect, we have opted for the two-step procedure.

Appendix D: Subluminal versus superluminal scenario in Simulation B and C

Fig. D.1. Same as Fig. 3 but for Simularion B. See the caption of Fig. 3 for details.

Fig. D.2. Same as Fig. 3 but for Simularion C. See the caption of Fig. 3 for details.

Appendix E: Simulation A, B and C fluxes

Fig. E.1. Normalized observer frame synchrotron flux of Simulation A (blue circles), B (orange triangles), and C (green squares) as a function of frequency.

Appendix F: Simulation A emissivities

Fig. F.1. Normalized observer frame synchrotron emissivity in the radio band of a 10-code units slab centered on xz plane (Simulation A).

Fig. F.2. Same as Fig. F.1 but in the optical band.

Fig. F.3. Same as Fig. F.1 but in the X-ray band.

Appendix G: Simulation C emissivities

Fig. G.1. Normalized observer frame synchrotron emissivity in the radio band of a 10-code units slab centered on xz plane (Simulation C).

Fig. G.2. Same as Fig. G.1 but in the optical band.

Fig. G.3. Same as Fig. G.1 but in the X-ray band.

Appendix H: Effect of the viewing angle on the polarization degree

Fig. H.1. Polarization degree of Simulation A as a function of frequency for different viewing angles.

All Tables

All Figures

	Fig. 1. Map of the pressure (left) and Lorentz factor (right) in the xz plane (Simulation A). The main structures are indicated. The cells tagged as shock regions in the xz plane are highlighted in light blue.
In the text

	Fig. 2. Map of the Lorentz factor in the xz plane with the projection of the trajectory of 3 macroparticles (number 799, 3983, 7789, Simulation C). The trajectory color represents the macroparticle status: white when inside a shock region and black when outside. The cells tagged as shock regions in the xz plane are highlighted in light blue.
In the text

	Fig. 3. Left and middle panels: Distribution of the upstream magnetic inclination and shock speed in the upstream rest frame for Simulation A. Right panel: Comparison of these two quantities for each configuration. The dot color indicates the z coordinate of the macroparticle exiting point from the shock. The configurations inside the red region are superluminal.
In the text

	Fig. 4. Map of the Lorentz factor in the xz plane with the projection of the trajectory of a macroparticle (number 1135, Simulation B). The color of the trajectory represents the time step, see Fig. 5.
In the text

	Fig. 5. Time evolution of the normalized spectrum of the macroparticle in Fig. 4.
In the text

	Fig. 6. Normalized observer frame synchrotron emissivity in the radio band of a 10-code units slab centered on xz plane (Simulation B).
In the text

	Fig. 7. Same as Fig. 6 but in the optical band.
In the text

	Fig. 8. Same as Fig. 6 but in the X-ray band.
In the text

	Fig. 9. Polarization degree (top panel) and angle (bottom panel) of Simulation A (blue circles), B (orange triangles), and C (green squares). The polarization angles are measured clockwise from the projection of the jet axis onto the sky plane (solid gray line).
In the text

	Fig. 10. Map of the poloidal (left panel) and toroidal (right panel) component of the fluid rest frame magnetic field in the xz plane (Simulation B).
In the text

	Fig. A.1. Profiles of l(θ) and h(θ)/Γ.
In the text

	Fig. D.1. Same as Fig. 3 but for Simularion B. See the caption of Fig. 3 for details.
In the text

	Fig. D.2. Same as Fig. 3 but for Simularion C. See the caption of Fig. 3 for details.
In the text

	Fig. E.1. Normalized observer frame synchrotron flux of Simulation A (blue circles), B (orange triangles), and C (green squares) as a function of frequency.
In the text

	Fig. F.1. Normalized observer frame synchrotron emissivity in the radio band of a 10-code units slab centered on xz plane (Simulation A).
In the text

	Fig. F.2. Same as Fig. F.1 but in the optical band.
In the text

	Fig. F.3. Same as Fig. F.1 but in the X-ray band.
In the text

	Fig. G.1. Normalized observer frame synchrotron emissivity in the radio band of a 10-code units slab centered on xz plane (Simulation C).
In the text

	Fig. G.2. Same as Fig. G.1 but in the optical band.
In the text

	Fig. G.3. Same as Fig. G.1 but in the X-ray band.
In the text

	Fig. H.1. Polarization degree of Simulation A as a function of frequency for different viewing angles.
In the text