# Shock corrugation to the rescue of the internal shock model in microquasars
The single-scale MHD view

Patryk Pjanka Nordita, KTH Royal Institute of Technology and Stockholm
University, Hannes Alfvéns väg 12, SE-106 91 Stockholm, Sweden  Camilia
Demidem JILA, University of Colorado and National Institute of Standards and
Technology,
440 UCB, Boulder, CO 80309-0440, USA  Nordita, KTH Royal Institute of
Technology and Stockholm University, Hannes Alfvéns väg 12, SE-106 91
Stockholm, Sweden  Alexandra Veledina Department of Physics and Astronomy,
FI-20014 University of Turku, Finland  Nordita, KTH Royal Institute of
Technology and Stockholm University, Hannes Alfvéns väg 12, SE-106 91
Stockholm, Sweden

###### Abstract

Questions regarding energy dissipation in astrophysical jets are open to date,
despite of numerous attempts to limit the diversity of models. Some of the
most popular models assume that energy is transferred to particles via
internal shocks, which develop as a consequence of non-uniform velocity of the
jet matter. In this context, we study the structure and energy deposition of
colliding plasma shells, focusing our attention on the case of initially
inhomogeneous shells. This leads to formation of distorted (corrugated) shock
fronts – a setup that has recently been shown to revive particle acceleration
in relativistic magnetized perpendicular shocks. Our studies show that the
radiative power of the far downstream of non-relativistic magnetized
perpendicular shocks is moderately enhanced with respect to the flat shock
cases. Based on the decay rate of downstream magnetic field, we make
predictions for multiwavelength polarization properties.

Jets (870) – Shocks (2086) – Magnetohydrodynamical simulations (1966) –
Particle astrophysics (96) – X-ray binary stars (1811) – Active galactic
nuclei (16)

##  1 Introduction

Fast relativistic outflows – jets – are common to astrophysical environments
harboring accretion and / or explosions. Sources standing orders of magnitude
apart in energy dissipation, such as young stellar objects, accreting black
hole X-ray binaries, gamma-ray bursts, active galaxies and blazars, have all
been witnessed to launch such outflows (e.g., Goodson et al., 1997; Mirabel &
Rodríguez, 1999; MacFadyen & Woosley, 1999). The mechanisms behind jet
acceleration and fueling are, however, debated (Komissarov, 2011; Davis &
Tchekhovskoy, 2020). Jets are known to produce bright non-thermal emission,
often observed as power-law synchrotron and/or inverse Compton continua (e.g.,
review by Matthews et al., 2020), indicating continuously-operating energy
dissipation and particle acceleration mechanisms.

Jet kinetic energy can be converted to radiation via magnetic reconnection
events (Lyutikov, 2003), pair processes (Derishev et al., 2003; Stern &
Poutanen, 2008) or shocks in the interior of the jets (Pe’er, 2014). In the
latter scenario, the jet is thought to be ejected as inhomogeneous matter,
both in density and in velocity. This gives rise to discrete pockets of plasma
(shells) travelling through the jet at different speeds. Once a faster shell
catches up with a (previously ejected) slower one, they collide, and two
shocks are formed (under the right conditions, see Pe’er 2014): a forward
shock driven into the slower shell, and a reverse shock driven into the faster
one. These shocks could enable particle acceleration and, consequently, power
the non-thermal emission of the jet. After the collision, the shells are
typically assumed to merge inelastically and continue propagating as a single
shell (e.g., Beloborodov 2000; Malzac 2014, although see Kino et al. 2004).
They then encounter other shells and collide again, ensuring continuous
particle re-acceleration throughout the jet. This “internal shock model” has
been originally proposed in the context of the resolved inhomogeneous
structures, the so-called knots, in the jet of M87 galaxy (Rees, 1978). Later,
the model was generalised to jets observed from other astrophysical sources:
radio-loud quasars (Spada et al., 2001), gamma-ray bursts (Rees & Meszaros,
1994; Medvedev & Loeb, 1999; Beloborodov, 2000), and microquasars (Kaiser et
al., 2000; Malzac, 2014).

Such bright knots have been observed in a number of sources, supporting
localised energy dissipation scenarios (Matthews et al., 2020). The observed
collision of these knots in the radio galaxy 3C 264 (Meyer et al., 2015) gave
observational confirmation to the internal shock scenario in context of AGN
jets. In microquasars, internal shock models have proved to be successful in
explaining both spectral and timing properties (Vincentelli et al., 2019;
Péault et al., 2019). However, such modelling relies on a number of
assumptions about energy liberation and particle acceleration in shell
collisions, which need to be verified.

Different aspects of the flow structure and jet observables in the presence of
internal shocks have been considered: the geometrical and hydrodynamic
evolution (Kino et al., 2004; Mimica et al., 2005; Granot, 2012; Marino et
al., 2020; Rudolph et al., 2020), thermodynamics (Kobayashi & Sari, 2001;
Mimica et al., 2004; Graff et al., 2008; Pe’er et al., 2017), the impact of
magnetic fields (Fan et al., 2004; Mimica et al., 2007; Mimica & Aloy, 2010,
2012; Rueda-Becerril et al., 2014, 2015; Deng et al., 2015), radiative
processes (Graff et al., 2008; Böttcher & Dermer, 2010; Joshi & Böttcher,
2011; Jamil & Böttcher, 2012), and timing properties (Spada et al., 2001;
Mimica et al., 2005; Böttcher & Dermer, 2010; Jamil et al., 2010; Malzac,
2013, 2014; Drappeau et al., 2015), among others. While the internal shock
scenario is quite successful in representing a number of sources (e.g., Daigne
& Mochkovitch, 1998; Spada et al., 2001; Drappeau et al., 2015; Malzac et al.,
2018; Péault et al., 2019; Bassi et al., 2020; Marino et al., 2020), its main
problem (and its main criticism) lies in its apparent inefficiency of energy
conversion from kinetic to that of non-thermal particles (though this may
depend on model assumptions, as pointed out by Pe’er et al. 2017). It has been
shown (Begelman & Kirk, 1990; Sironi & Spitkovsky, 2009, 2011) that electrons
can consistently achieve non-thermal (power-law) distributions only in
parallel shocks, where the magnetic field is nearly parallel to the shock
normal. Once the obliquity of magnetic field increases, particle acceleration
is inhibited.

Recently, a potential solution to this issue has been proposed. Efficient
particle acceleration can be recovered for the case of shocks corrugated
(rippled) by upstream inhomogeneities (Demidem et al., 2022). Such
inhomogeneities may be produced by inhomogeneous magnetic field configurations
(e.g., Cerutti & Giacinti, 2020), or as a result of operating kinetic
turbulence (e.g., Zhdankin et al., 2017). Interestingly, they were also seen
to self-consistently arise in kinetic electron-ion simulations (Ligorini et
al., 2021). These findings can recover internal shocks as an efficient
mechanism to produce nonthermal particles in relativistic outflows.

The structure behind (downstream of) shocks corrugated by upstream/downstream
perturbations is much more complex than for planar shocks in homogeneous
media. Demidem et al. (2018) described these perturbations in terms of
eigenmodes of linear MHD (i.e., entropy, Alfvén, and magnetosonic waves) and
studied the response of a fast perpendicular shock to each mode, one at a
time. Their MHD simulations illustrate how an incoming upstream plane wave
induces shock corrugation and is transmitted into a collection of modes in the
downstream (see, e.g., their first figure or, for corrugation induced by the
reflection of downstream perturbations, Fig. 6 of Lyutikov et al. 2012).
Demidem et al. (2022) used a kinetic approach (enabling study of particle
energization from first principles) and investigated possible implications for
particle acceleration in relativistic magnetized perpendicular shocks. While
the process is inhibited for shocks in initially homogeneous media, they find
that shock rippling and downstream turbulence induced by (upstream) harmonic
density perturbations can lead to generation of energetic particles.

Inspired by these results, the present work proposes to analyze the properties
of a collision of two-dimensional (2D) corrugated shells and how its
observational signatures may differ from the idealized (one-dimensional, or,
1D) case of flat shells. We perform relativistic magnetohydrodynamic (MHD)
simulations of colliding shells, corrugated by ambient density or pressure
inhomogeneities, and follow the evolution of density, magnetic field, and
energy dissipation (as probed by synchrotron emission). We find that the
quantities are moderately affected by the 2D-structure of the shocks, relative
to the 1D case with flat colliding shells. Hence, we conclude that the
previous predictions of analytical models and hydrodynamic simulations hold
for the case of acceleration-efficient corrugated shocks.

##  2 Methods

The one-dimensional behaviour of colliding shells, including the effects of
shell magnetization, initial internal energy, and other parameters, have been
investigated in exquisite detail by previous works (e.g., Kobayashi & Sari,
2001; Fan et al., 2004; Kino et al., 2004; Mimica et al., 2004, 2005, 2007;
Jamil et al., 2010; Mimica & Aloy, 2010; Joshi & Böttcher, 2011; Mimica &
Aloy, 2012; Pe’er et al., 2017; Rudolph et al., 2020). Thus, here we treat the
one-dimensional behavior of the shells as well established, and instead focus
on the difference between a 1D shell collision and the matching (see Sect.
2.4) 2D corrugated case.

In each (1D/2D) case, we simulate a collision of two identical dense shells,
which are magnetized uniformly, with magnetic field lines being perpendicular
to their motion (corresponding to the so-called perpendicular shock
configuration), and moving through a low-density, unmagnetized ambient medium.
We simulate the encounter in the zero-momentum frame of the system. For each
test case, we run two simulations: a one-dimensional one (henceforth 1D) with
no shell structure perpendicular to the shell motion (Sect. 2.3), and a two-
dimensional simulation (2D) with shell surface corrugation (Sect. 2.4).

###  2.1 Numerical setup

We use the grid-based astrophysical MHD code Athena 4.2111Athena 4.2 is
publicly available at https://github.com/PrincetonUniversity/Athena-Cversion.
A fork of this repository augmented to run and process the simulations
presented in this paper is available at
https://github.com/ppjanka/Corrugated_internal_shocks. (Stone et al., 2008), a
Godunov code using constrained transport (CT) method to enforce the zero-
divergence condition on magnetic fields.

Our simulations are performed in special relativistic MHD (Beckwith & Stone,
2011), in Cartesian coordinates. The boundary conditions in the direction
parallel to the shell motion (along the jet) are set as “free outflow”
(variables at the domain edge are copied to the ghost cells), while those
perpendicular to shell motion are periodic. The equation of state is adiabatic
with the adiabatic index of $\gamma_{\rm ad}=4/3$. We use 2nd-order
reconstruction in primitive variables, the van Leer unsplit time integrator
(Stone & Gardiner, 2009), and HLLC (Toro, 1999; Mignone et al., 2009) as our
Riemann solver, with $\textrm{CFL}=0.025$.

For our grid resolution, we require that i) for each corrugation method (see
Sect. 2.4) we resolve the corrugation width $w$ of the shock by at least $8$
cells (parallel to the mean shock normal, see Sect. 2.5 and Fig. 3) and ii)
that the resolution is nearly identical in both grid directions (i.e., our
cells are nearly square-shaped). This results in the simulation resolution of
$8192\times 1024$ for a box of $40$ light seconds $\times$ $6$ light seconds
(hereafter, lt-sec). For 1D simulations, we use the same resolution along the
jet (for fair comparison) and collapse the perpendicular direction to a single
cell.

###  2.2 Choice of parameter space

Microphysics
---
$\gamma_{\rm ad}$ | adiabatic index | $4/3$
Geometry
Box size (lt-sec) | $40\times 6$
1D box size (cells) | $8192\times 1$
2D box size (cells) | $8192\times 1024$
Distance between shells (lt-sec) | $5$
Shell width (lt-sec) | $1$
Shell fluid
$\rho_{\rm sh}$ | Density | $10^{-13}\textrm{ g}/\textrm{cm}^{3}$
$\sigma$ | Magnetization | $0.1$
$\beta$ | Plasma $\beta$ | $1$
$v_{\rm sh}$ | Bulk speed (collision frame) | $0.1c$
Ambient medium
$\rho_{\rm amb}$ | Density | $10^{-3}\rho_{\rm sh}$
$P_{\rm amb}$ | Pressure | $10^{-17}(\textrm{g}/\textrm{cm}^{3})c^{2}$
$B_{\rm amb}$ | Magnetic field strength | 0
Synchrotron flux
$D$ | Distance |  $8$ kpc
$R$ | Size of the emitting region |  $6\times 10^{8}$ cm
| (geometrical depth) |
$\Gamma_{\rm j}$ | Jet bulk Lorentz factor | $2$
$i$ | Jet axis inclination | $60\degree$
|  | or $30\degree\sim 1/\Gamma_{\rm j}$
Table 1: Parameters of our model. Where applicable, initial values are given –
note that these may be altered by the corrugation procedure (both for 1D and
2D runs; as explained in Sect. 2.4). For synchrotron flux calculations, we use
microquasar-relevant values as our fiducial model; we also discuss the
application of our results to blazars in text. $c$ denotes the speed of light.

To make our calculations more readable, we set our problem parameters and
units to match the environment of an existing source (at least to an order of
magnitude) and list them in Table 1. While our calculations are relevant for
any relativistic jet where emission is synchrotron-dominated, we provide
numerical values for the microquasar GX 339$-$4\. We briefly discuss a
comparison with blazars, another prime candidate for application of the
internal shock model, at the end of this section.

GX 339$-$4 is a well-observed low-mass black hole X-ray binary, with a short
duty cycle of X-ray outbursts of $\sim 2-3$ years (see, e.g., Corbel et al.,
2013; Connors et al., 2019). While its orbital period can be measured
precisely at $1.76$ days (Hynes et al., 2003; Levine & Corbet, 2006; Heida et
al., 2017), the distance, compact object mass (and nature), and inclination
are difficult to specify from observations. Heida et al. (2017) constrain
these to the ranges of $D\gtrsim 5$ kpc, $M\in(2.3,9.5)M_{\odot}$, and
$i\in(37\degree,78\degree)$, respectively. For the purpose of this work, we
will adopt the values of $D\sim 8$ kpc and $M\sim 10$ $M_{\odot}$. We also
adopt a jet opening angle of $\theta_{\rm j}\sim 2.3\degree$ and jet bulk
Lorentz factor of $\Gamma_{\rm j}\sim 2$ (as assumed by Malzac et al. 2018, a
value typical of microquasar jets, e.g., Tetarenko et al. 2019; Zdziarski et
al. 2022b). We performed our calculations for two different inclinations:
$i=30\degree\sim 1/\Gamma_{\rm j}$ (a head-on jet, corresponding to a blazar-
like scenario) and $i=60\degree$ (a jet seen from its side, more typical of
X-ray binary jets). However, we find that the relative differences between
corrugated and non-corrugated runs (on which we base our analysis, see Sect.
2.5) are very weakly dependent on inclination. Thus, we only report the
results for $i=60\degree$ in this manuscript.

If we scale our simulations to microquasar conditions, our region of interest
would be radiating in the IR. Based on Malzac et al. (2018), this corresponds
to $z\sim 10^{3}-10^{6}r_{\rm g}$ (see their fig. 4), where $r_{\rm
g}=1.5\times 10^{6}$ cm for a BH of $10M_{\odot}$. Shidatsu et al. (2011)
estimated the jet IR-emitting region size at $R\sim 6\times 10^{8}$ cm.
Assuming $i=2.3\degree$, this gives $z\sim 1.5\times 10^{10}\textrm{ cm}\sim
10^{4}r_{\rm g}$, which we adopt here.

In continuous jet models, applicable up to the radio-emitting regions far in
the jet, it is often assumed that the high-energy particle and magnetic field
energy densities are in equipartition (e.g., Hirotani 2005, Zdziarski et al.
2012, Drappeau et al. 2015; see, however, Zdziarski et al. 2015). For
simplicity, we also follow this assumption here. To obtain the remaining
plasma parameters, we use existing theoretical models of the binary from the
literature to provide lower and upper constraints. We show detailed
derivations and the resulting ranges in Appendix A. Informed by these results,
we adopt the rest-frame fluid density of $\rho\sim 10^{-13}\textrm{
g}/\textrm{cm}^{3}$ and magnetization of $\sigma\sim 0.1$.

As noted above, while we report numerical values for a microquasar, our models
can also be interpreted in terms of blazar sources. This would require a
number of changes in Table 1, most of which would result in similar
enhancements or decreases of synchrotron emission for both the 1D and 2D cases
(i.e., change of normalization). We are, however, only interested in the
difference between the paired 1D/2D models. We find this difference to be
mainly influenced by two parameters. Magnetization $\sigma$ contains
information about how dominant magnetic fields are in the dynamics of the
problem. For large $\sigma$, magnetic fields are able to prevent corrugation,
leading to more similar 1D and 2D results. Fortunately, magnetization of the
ballistic parts of a blazar jet can be estimated to $\sigma_{\infty}\sim
0.1-0.01$ (see Tchekhovskoy et al., 2009; Zdziarski et al., 2015, and
references therein), so conditions in blazar jets are similar to the ones
presented here in that respect. The second parameter of relevance is the bulk
Lorentz factor of the jet $\Gamma_{\rm j}$. Faster jets likely contain a
broader distribution of shell speeds, with high-speed collisions (large
$v_{\rm sh}$) more likely. In the models presented in our work, collisions are
driven by adiabatic expansion of the shells into the ambient medium between
them. At faster relative speeds, this expansion may become less relevant to
the collision outcome. If shell corrugation is caused by pressure imbalance
within the shell (see our pressure-corrugated runs, Sect. 2.4), this will
result in flatter shocks, and more similar 1D and 2D runs (i.e., less effect
of corrugation in the synchrotron emission).

###  2.3 1D initial conditions

Figure 1: A diagram of our initial conditions and numerical setup (not to
scale) before corrugation is applied.

Throughout the paper, the fluid quantities (density $\rho$ and thermal
pressure $P$) are given in the fluid rest-frame, while the velocity $v$,
Lorentz factor $\Gamma$, and magnetic field $B$ are evaluated in the
simulation frame (zero-momentum frame of the shell collision). We denote the
coordinate along the jet, parallel to the movement of the shells, as $x$ and
the perpendicular one as $y$.

The initial conditions for non-corrugated (1D) runs are depicted in Fig. 1.
The domain is composed of 5 regions delimited by the shell edges: $x_{\rm
L1}$, $x_{\rm L2}$ for the left shell and $x_{\rm R1}$, $x_{\rm R2}$ for the
right shell, with the subscript “1” closer to the center of the domain and “2”
further away. Note that our setup is symmetric with respect to the plane of
the (eventual) shell collision ($x=0$, a vertical line at the center of our
domain). Within the shells, density is $\rho_{\rm sh}$, pressure $P_{\rm sh}$,
velocity $v_{\rm sh}$ is directed towards the center of the domain, and the
magnetic field is set to $B_{\rm sh}\vec{e}_{y}$. Outside the shells, density
equals $\rho_{\rm amb}$, pressure $P_{\rm amb}$ and magnetic field is 0. The
velocity field in the ambient medium is initialized according to:

$\vec{v}=\left\\{\begin{array}[]{cc}v_{\rm sh}\vec{e}_{x},&x\leq x_{\rm
L1},\\\ 0,&x\in(x_{\rm L1},x_{\rm R1}),\\\ -v_{\rm sh}\vec{e}_{x},&x\geq
x_{\rm L1}.\\\ \end{array}\right.$ (1)

Within the shells, the magnetic field strength $B_{\rm sh}$ and pressure
$P_{\rm sh}$ are set through the plasma $\beta$ and $\sigma_{B}$ parameters:

$B_{\rm sh}=\sqrt{8\pi\rho_{\rm sh}\sigma_{B\rm,sh}}\Gamma_{\rm sh}c,$ (2)
$P_{\rm sh}=\beta\frac{B_{\rm sh}^{2}}{8\pi\Gamma_{\rm sh}^{2}}.$ (3)

Numerically, we set our initial conditions to $\rho_{\rm sh}=10^{-13}\textrm{
g}/\textrm{cm}^{3}$, $\beta=1$, $\sigma=0.1$ (cf., Sect. 2.2 and Table 1),
$v_{\rm sh}=0.1c$ (where $c$ denotes the speed of light).

###  2.4 Shock corrugation methods

To produce shock corrugation, we modulate values of parameters from the
previous paragraph in the direction perpendicular to shell motion. In studies
of corrugated shocks, various (typically dynamical) methods are used to induce
corrugation, or bending, of the shock front (see, e.g., Demidem et al., 2018).
We investigate two possibilities: density and pressure modulation. In real
astrophysical sources, such shock-rippling perturbations can follow very
complex spatial dependencies. However, they can always be decomposed into
harmonic components. Thus, as a first step, we use sinusoidal modulation to
represent single harmonics of such shock-rippling density or pressure
perturbations in the fluid, in a fashion similar to that of Demidem et al.
(2018, 2022).

An important aspect of each implementation lies in ensuring that each 1D/2D
pair of simulations (with matching global parameters) is only distinguished by
the 2D structure of the latter, and not by, e.g., the initial total
synchrotron flux. As a result, each modulation in the initial conditions
creating a 2D (corrugated) simulation may also include a change in the
parameters of the corresponding 1D case (i.e., modify its parameters:
$\rho_{\rm amb}$, $\rho_{\rm sh}$, $P_{\rm amb}$, etc.). This is done to
ensure that the initial lab-frame total energy (integrated over the simulation
domain),

$\epsilon_{\rm tot}=\int\left[\Gamma^{2}\left(\rho c^{2}+\frac{\gamma_{\rm
ad}}{\gamma_{\rm ad}-1}P\right)-P-\Gamma\rho c^{2}\right]dV,$ (4)

shell momentum,

$p_{\rm sh,tot}=\int_{\rm sh}\left[\Gamma^{2}\left(\rho
c^{2}+\frac{\gamma_{\rm ad}}{\gamma_{\rm
ad}-1}P\right)\frac{v}{c^{2}}\right]dV,$ (5)

and synchrotron power (see Sect. 2.5),

$F_{\rm syn,tot}\mathchoice{\mathrel{\vbox{
\offinterlineskip\halign{\hfil$#$\cr\displaystyle\propto\cr\kern
2.0pt\cr\displaystyle\sim\cr\kern-2.0pt\cr}}}}{\mathrel{\vbox{
\offinterlineskip\halign{\hfil$#$\cr\textstyle\propto\cr\kern
2.0pt\cr\textstyle\sim\cr\kern-2.0pt\cr}}}}{\mathrel{\vbox{
\offinterlineskip\halign{\hfil$#$\cr\scriptstyle\propto\cr\kern
2.0pt\cr\scriptstyle\sim\cr\kern-2.0pt\cr}}}}{\mathrel{\vbox{
\offinterlineskip\halign{\hfil$#$\cr\scriptscriptstyle\propto\cr\kern
2.0pt\cr\scriptscriptstyle\sim\cr\kern-2.0pt\cr}}}}B^{\frac{p+5}{2}},$ (6)

of paired 1D/2D runs are identical (and thus, the runs are comparable). In the
equations above, $\Gamma$ is the local fluid Lorentz factor (in the simulation
frame), $\rho$ – the local fluid-frame mass density, $P$ – fluid-frame thermal
pressure, and $B$ is the simulation-frame magnetic field strength.
Conservation of quantities (4)-(6) can be achieved by ensuring that the
averages of $\rho$, $P$, and $B$ over the $y$ coordinate match between
corresponding 1D/2D cases.

####  2.4.1 Density modulation

In our first method of inducing shell corrugation, that we will refer to as
“density corrugation”, we change the density of ambient medium between shells
to exhibit “ridges” (sinusoidal density modulation oblique with respect to the
shell surface) that bend the surface of each shell as it sweeps their mass. To
ensure that the paired 1D run has initially identical global properties
(4)-(6), we also increase the ambient medium density $\rho_{\rm amb}$ to match
the mean density in the “ridged” zone in 2D run. In other words, for the 2D
case we set:

$\rho_{\rm 2D}=\left\\{\begin{array}[]{ll}\xi\rho_{\rm amb},&x<x_{\rm L2},\\\
\rho_{\rm sh},&x\in[x_{\rm L2},x_{\rm L1}],\\\ \begin{array}[]{rl}\rho_{\rm
amb}+\frac{1}{2}&(A_{\rho}\rho_{\rm sh}-\rho_{\rm amb})\\\
&\times\mathcal{M}_{\rho}(x,y)\end{array},&x\in]x_{\rm L1},x_{\rm R1}[,\\\
\rho_{\rm sh},&x\in[x_{\rm R1},x_{\rm R2}],\\\ \xi\rho_{\rm amb},&x>x_{\rm
R2},\end{array}\right.$ (7)

where $A_{\rho}$ is the corrugation amplitude,
$\displaystyle\xi=\frac{1}{2}\left(1+A_{\rho}\frac{\rho_{\rm sh}}{\rho_{\rm
amb}}\right)$, and

$\mathcal{M}_{\rho}(x,y)=\cos\left(2\pi(n_{x}\frac{x-x_{\rm L1}}{x_{\rm
R1}-x_{\rm L1}}\right.\\\ \left.+n_{y}\frac{y-y_{\rm min}}{y_{\rm max}-y_{\rm
min}})\right)+1,$ (8)

with $(n_{x},n_{y})=(1,2)$. For the 1D case (to conserve quantities in Eqs. 4
and 5), we set

$\rho_{\rm 1D}=\left\\{\begin{array}[]{ll}\xi\rho_{\rm amb},&x<x_{\rm L2},\\\
\rho_{\rm sh},&x\in[x_{\rm L2},x_{\rm L1}],\\\ \xi\rho_{\rm amb},&x\in]x_{\rm
L1},x_{\rm R1}[,\\\ \rho_{\rm sh},&x\in[x_{\rm R1},x_{\rm R2}],\\\
\xi\rho_{\rm amb},&x>x_{\rm R2}.\end{array}\right.$ (9)

Note that we only adjust the ambient medium conditions, and leave the fluid
within the shells unchanged. Thus, 1D/2D pairs of runs with different
$A_{\rho}$ have different density contrast ratios. While the runs with
different $A_{\rho}$ are in effect less easily comparable in absolute terms
(and so we only use relative values in all our comparisons), it allows us to
focus on the effect of corrugation.

Overall, we perform density-corrugation simulations (1D/2D pairs) for
$A_{\rho}$ values of $0.01$, $0.02$, $0.05$, $0.10$, $0.20$, $0.50$, $0.75$,
and $1.00$.

####  2.4.2 Shell pressure modulation

In the second method of inducing shell corrugation, henceforth referred to as
“pressure corrugation”, we implement a pressure modulation within the shell.

$P_{\rm sh,2D}=P_{\rm sh}+\frac{1}{2}A_{P}(P_{\rm sh}-P_{\rm amb})\\\
\times\cos\left[2\pi\left(n_{x}\frac{x-x_{\rm S2}}{x_{\rm S1}-x_{\rm
S2}}+n_{y}\frac{y-y_{\rm min}}{y_{\rm max}-y_{\rm min}}\right)\right],$ (10)

where $S\in\\{L,R\\}$ is the shell index, $(n_{x},n_{y})=(1,2)$, and $A_{P}$
is the “corrugation amplitude”. As the conditions in Eq. (10) do not change
our original mean internal energy, shell momentum, or synchrotron flux, other
initial conditions remain unchanged.

We perform pressure-corrugation simulations (1D/2D pairs) for $A_{P}$ values
of $0.10$, $0.20$, $0.50$, $0.75$, and $1.00$.

###  2.5 Diagnostics

####  2.5.1 Measuring corrugation

Figure 2: An example measurement of shock corrugation aspect ratio ($a_{\rm
sh}=w/\lambda$), for a pressure corrugated case with $A_{P}=0.2$. Top: a
density contour plot with shock detection marked as red lines and aspect ratio
measurement depicted with the white curve and annotations (see text). Bottom:
shock width $w$ as a function of time (with respect to collision time $t_{\rm
coll}$). The dashed black and red lines show the moments of collision and
maximal corrugation, respectively. The final shock width $w_{\rm max}$ is
measured at the red dashed line. Figure 3: Scaling between the shock aspect
ratio and the initial conditions’ corrugation parameter for density
corrugation (blue) and pressure corrugation (orange).

In our models, shell corrugation is induced in the initial conditions, via the
modulation fractions $A_{\rho}$ and $A_{P}$. These parameters tell us little
about the eventual shock geometry. Therefore, we explicitly measure the shape
of the corrugated shock we obtain in each simulation. As our diagnostic, we
use the shock’s “aspect ratio” (see Fig. 2), the ratio of its width $w$, the
distance between maximum corrugation down- and upstream, to its wavelength
$\lambda$, set by the parameters $(n_{x},n_{y})=(1,2)$ in Eqs. (8–10).

From the 1D version of a given run, we extract the maximal density for each
frame, and note the time of maximal compression – which we consider the
“collision time” $t_{\rm coll}$. For the 2D version, we inspect the frames
between $t_{\rm coll}-t_{\rm cross}$ and $t_{\rm coll}+nt_{\rm cross}$ (with
$n=2,4$ for density and pressure corrugation, respectively), where $t_{\rm
cross}$ is the light-crossing time of the initial shell width. For each or
those frames, we measure the instantaneous shock width:

  1. 1.

In each frame, the center of the frame is surrounded from both sides by two
shocks, either driven by shells propagating through the ambient medium
($t<t_{\rm coll}$), or driven into each shell by the collision ($t>t_{\rm
coll}$). We localize these shocks using the velocity gradient. $x$-coordinates
of points with $|\partial v_{x}/\partial x|$ larger than a threshold value are
averaged along the $x$ direction separately for the left and right halves of
the simulated domain, giving us the precise locations of both shocks at each
$y$ (shown by red curves in Fig. 2).

  2. 2.

At each $y$, we average the $x$-position of these two shocks horizontally, to
obtain the white stripe (along the center) in Fig. 2.

  3. 3.

We then collect the $x$ coordinates of all points within the white stripe (one
per each $y$ value) and produce their histogram. To discard the outliers, we
take the distance between the maxima of this histogram to be the width of the
shock in a given frame.

In all our simulations, we find that the shock width measured this way
decreases before the collision, reaches a minimum at $t_{\rm coll}$, then
increases and reaches a maximum, after which the width continuously decreases
as the post-shock regions evolve. We take the post-collision maximal shock
width $w_{\rm max}$ as the shock width for the given 1D/2D configuration
(i.e., for the given $A_{\rho}$ or $A_{P}$). The corrugation wavelength is
always $\Delta y/2$ (where $\Delta y$ is the simulated box’s height; cf. Eqs.
(8) – (10)), giving us the aspect ratio of:

$a_{\rm sh}=2w_{\rm max}/\Delta y.$ (11)

We plot the relations between $a_{\rm sh}$ and $A_{\rho}$, $A_{P}$ in Fig. 3.

####  2.5.2 Synchrotron emission

The main deliverable of internal shock models lies in their accurate
prediction of jet emission which, in the leptonic model adopted here, is
powered by the non-thermal electron population. Of course, an MHD-only picture
does not include non-thermal particles, and thus a number of assumptions must
be adopted to extract the desired observables. We select our assumptions in
line with typical prescriptions from the literature (cf., e.g., Spada et al.,
2001; Graff et al., 2008; Böttcher & Dermer, 2010; Jamil et al., 2010; Malzac,
2014; Rueda-Becerril et al., 2014; Malzac et al., 2018). At each location
within the jet, the local non-thermal electron energy distribution is assumed
to be a power law $dN/d\gamma\propto\gamma^{-p}$ with a slope of $p=2.5$,
extending between the Lorentz factors $\gamma_{\rm min}=10$ and $\gamma_{\rm
max}=10^{6}$. Non-thermal electrons are taken to be in energy equipartition
with the magnetic field (i.e., with equipartition parameter
$\zeta_{e}=U_{e}/U_{B}=1$, where $U_{e}=\int(\gamma
m_{e}c^{2}\frac{dN}{d\gamma})d\gamma$ and $U_{B}=B^{2}/(8\pi)$).

We stress that, as a consequence of the assumptions above, the impact of
thermal particles on the acceleration and radiative processes is assumed to be
negligible – thermal particles (represented by the mass density in MHD
simulations) only affect the fluid dynamics of the model. Particle
acceleration is assumed to proceed until there is equipartition between the
non-thermal particles and the magnetic field. Thus, the number density of non-
thermal particles is a direct function of the magnetic field energy density.
This is regardless of the MHD mass density, which is assumed to be
sufficiently large to provide seed particles for the non-thermal distribution
without itself being significantly depleted. In turn, we also assume that
synchrotron emission and self-absorbtion are both completely dominated by non-
thermal particles – we caution that this last assumption may not be applicable
in some blazars, where self-absorbtion by thermal electrons may be important.

For our results, we report the synchrotron flux per unit jet surface area in
the observer’s frame (e.g., Rybicki & Lightman, 1979):

$\frac{dF_{\nu}}{dS}=\frac{\delta^{2}\Gamma_{\textrm{fl}\rightarrow\textrm{obs}}}{2D^{2}}f_{v}\bar{I}_{\bar{\nu}},$
(12)

where
$\delta=[\Gamma_{\textrm{fl}\rightarrow\textrm{obs}}(1-\vec{v}_{\textrm{fl}\rightarrow\textrm{obs}}\cdot\vec{n}/c)]^{-1}$
is the Doppler factor, $\vec{n}$ is the direction towards the observer,
$\Gamma_{\textrm{fl}\rightarrow\textrm{obs}}=1/\sqrt{1-v^{2}_{\textrm{fl}\rightarrow\textrm{obs}}/c^{2}}$,
$f_{\nu}=1$ is the volume filling factor of non-thermal electrons,
$\nu=\bar{\nu}\delta$ is the observer frame photon frequency linked to the
jet-frame frequency $\bar{\nu}$, and $\bar{I}_{\bar{\nu}}$ is the jet-frame
synchrotron intensity. In order to determine the values of the fluid
quantities in the observer’s frame from their rest-frame values, one needs to
account for both the velocity $\vec{v}=(v_{x},v_{y})$ of the fluid in the
simulation frame (zero-momentum frame of the shell collision) and the speed of
the simulation frame in the observer frame $\vec{v_{\rm j}}=v_{\rm
j}\vec{e}_{x}$ (bulk velocity of the jet). We do so by performing Lorentz
transformations of the fluid velocity:

$\vec{v}_{\textrm{fl}\rightarrow\textrm{obs}}=\frac{v_{x}+v_{\rm
j}}{1+v_{x}v_{\rm j}/c^{2}}\vec{e}_{x}+\frac{\sqrt{c^{2}-v_{\rm
j}^{2}}v_{y}}{c^{2}+v_{x}v_{\rm j}}\vec{e}_{y}.$ (13)

We consider the total IR-opt flux $dF_{\rm syn}/dS$, which we take to
correspond to the flux integrated over the frequency range $300$ GHz – $3$ PHz
range. For our assumed electron distribution, the local intensity can be
written as (e.g., Rybicki & Lightman, 1979; Beloborodov, 2000; Malzac, 2014):

$\bar{I}_{\bar{\nu}}=\frac{j_{\bar{\nu}}}{\alpha_{\bar{\nu}}}\left(1-e^{-\alpha_{\bar{\nu}}R}\right),$
(14)

where the synchrotron emission and extinction coefficients are:

$j_{\bar{\nu}}=K_{j}\zeta_{e}B^{\frac{p+5}{2}}\bar{\nu}^{-\frac{p-1}{2}},$
(15) $K_{j}=\frac{\sqrt{3}e^{3}i_{\gamma}}{16\pi^{2}m_{e}^{2}c^{4}(p+1)}\\\
\times\Gamma\left(\frac{3p+19}{12}\right)\Gamma\left(\frac{3p-1}{12}\right)\left(\frac{m_{e}c}{3e}\right)^{-\frac{p-1}{2}},$
(16)
$\alpha_{\bar{\nu}}=K_{\alpha}\zeta_{e}B^{\frac{p}{2}+3}\bar{\nu}^{-\frac{p+4}{2}}$
(17)
$K_{\alpha}=\frac{\sqrt{3}e^{3}i_{\gamma}}{64\pi^{2}m_{e}^{3}c^{4}}\left(\frac{3e}{2\pi
m_{e}c}\right)^{\frac{p}{2}}\\\
\times\Gamma\left(\frac{3p+2}{12}\right)\Gamma\left(\frac{3p+22}{12}\right),$
(18)

where $m_{e}$ and $e$ are the electron mass and unit charge, respectively,
and, for convenience, we define (following Malzac, 2014):

$i_{\gamma}^{-1}=\bar{\gamma}_{e}\int_{\gamma_{\rm min}}^{\gamma_{\rm
max}}\gamma^{-p}d\gamma$ (19)

with $\bar{\gamma}_{e}$ being the mean electron Lorentz factor. For the
optical depth $\alpha_{\bar{\nu}}R$, we select $R=6\times 10^{8}$ cm to match
the size of the IR-emitting core of the jet (Shidatsu et al., 2011). In other
words, we report our simulated box’s flux as if it were located along the jet
axis and observed with a line of sight at an angle of $1/\Gamma_{\rm j}$ to
the jet axis in the observer frame (i.e., perpendicular to the jet axis in the
comoving frame).

##  3 Results

A typical sequence of events during one of the 2D simulations is shown in Fig.
4. As the shells are over-pressurized with respect to their surrounding
ambient medium, the dynamics of the system is driven by adiabatic expansion of
the shells (for a thorough study of the dependence of shell collisions on the
initial internal energy see Pe’er et al. 2017). The shells start expanding
into the ambient medium with speeds in excess of their bulk velocities and
evolve under influence of the corrugation mechanism (see Sect. 2.4). Their
surface ripples in a linear manner, resulting in a near-sinusoidal shape. As
the shells collide, two shocks are driven into the shell material (one moving
left into the left shell and one moving right into the right one) and a thin
line of (unmagnetized) compressed ambient medium is trapped at a contact
discontinuity between the two post-shock regions. The post-shock material’s
density and magnetic field strength are both enhanced, with the latter causing
strong synchrotron emission (see Sect. 3.1 for further details). Note that,
due to our assumption of energy equipartition between the non-thermal
particles and the magnetic field, synchrotron emission is completely
determined by the latter. As the collision material continues to evolve, it
enters a non-linear phase, during which the central part of the domain becomes
disrupted. In the bottom two rows of Fig. 4, development of the Kelvin-
Helmholtz instability can be seen at $x\sim 0$ to mediate this process. While
the box-integrated total synchrotron flux remains high (albeit oscillating
with time), the emission becomes spread over the growing post-shock region,
and the local flux density dims considerably (as seen in the lower right
corner of Fig. 4).

Figure 4: A representative selection of snapshots of shell collision
evolution, taken from the density-corrugated (2D) simulation with
$A_{\rho}=0.2$ observed at the inclination of $i=60\degree$. Left column:
fluid-frame density, middle column: fluid-frame magnetic field strength, right
column: observer-frame synchrotron flux per unit emitting surface (log-scale).
Each row corresponds to a different time as annotated on the left (cf. Fig.
5). The linear and non-linear evolution of shell corrugation can be seen with
the progression of time. Note that, for the density corrugation case shown
here, only the inner shell edges (ones facing the other shell) are corrugated,
while the outer shell edges remain flat as they propagate outwards into the
ambient medium.

###  3.1 Corrugation-driven synchrotron flux enhancement

Figure 5: Top row: Light curves (mean synchrotron flux per surface emitting
area) for all of our simulations viewed at inclination of $i=60\degree$.
Middle and bottom rows: the difference between corrugated (2D) and non-
corrugated (1D) matched simulations shown as relative difference in
synchrotron flux integrated over the whole simulation domain (i.e., $2(F_{\rm
syn-tot,2D}-F_{\rm syn-tot,1D})/(F_{\rm syn-tot,2D}+F_{\rm syn-tot,1D})$). The
curves are color-coded by corrugation amplitudes $A_{\rho},A_{P}$ (see Fig. 3
for translation into shock aspect ratios).

We investigate how corrugation affects synchrotron flux from the shell
collision site. Fig. 5 shows the light curves integrated over our simulation
boxes (top row) and the relative difference between synchrotron fluxes of the
matched corrugated (2D) and non-corrugated (1D) cases. We notice two effects
introduced by shell corrugation.

There is a peak-then-dip pattern in the relative difference curves around
collision time ($\sim 3-8$ s). This is caused by the difference in collision
time. Due to corrugation, some regions of the 2D shell are travelling ahead of
the 1D shell position, and so the shells start colliding earlier in the
corrugated (2D) case – thus giving more emission than the 1D case just before
the 1D collision time (when the 2D collision happens), and less emission just
afterwards (which is at the 1D collision time, but already after the 2D
collision peak).

Figure 6: Effect of corrugation on the total synchrotron emission. Dependence
of the relative difference in the total synchrotron emission (flux integrated
over the simulated domain and simulation time) between the corrugated (2D) and
non-corrugated (1D) case, as a function of the shock aspect ratio (see Sect.
2.5 and Fig. 3). The point colors match the corresponding curves for each run
and corrugation method in Fig. 5. Note that, as a ratio of synchrotron fluxes,
this relation is independent of viewing inclination.

A more interesting result is visible in the long-term behavior of the relative
synchrotron flux differences. While the curves oscillate after collision time,
there is a clear overall enhancement of the synchrotron flux after shell
collision in the 2D case222Note that the outflow boundary conditions in the
direction along the flow remain causally disconnected from the center
throughout the simulation (and are therefore insensitive to the collision
geometry). Thus, this is a genuine effect of corrugation and not an artifact
caused by boundary conditions. – especially visible for the pressure-
corrugated case (bottom row of Fig. 5). To quantify this effect, we integrate
the light curves over simulation time, to obtain the total synchrotron
emission (per unit surface) for each simulation. We then subtract the 1D from
the 2D result and divide by their mean to obtain a total relative enhancement
of emission with respect to the 1D case. We plot these differences in Fig. 6.
We find that the differences can reach $14\%$ in our parameter range. They
also scale as a power law with the shock aspect ratio, with a slope of $1.5\pm
0.3$ and $2.38\pm 0.02$ for density- and pressure-corrugated runs,
respectively. Notably, the two curves for density- and pressure-corrugated
simulations match almost exactly. Therefore, the observed synchrotron
enhancement is independent of the corrugation mechanism and only sensitive to
the actual shock geometry.

Figure 7: Investigation of the cause of long-term synchrotron emission
enhancement. The curves show the relative synchrotron enhancement for the
corrugated (2D) simulation with respect to the matched uncorrugated one (1D;
cf. the bottom panels of Fig. 5) for the pressure corrugated case with
$A_{P}=0.5$. To calculate the emission in the corrugated case, the magnetic
field was either kept unchanged (blue curve), replaced with the magnetic field
from the uncorrugated (1D) case (orange), or the value was scaled down so that
the mean magnetic field strength matches that of the 1D simulation, with the
2D structure left unchanged (green). It can be seen that both the 2D structure
and an overall magnetic field enhancement contribute to the enhanced
synchrotron emission for the corrugated shell collisions.

Under the approximations adopted to calculate the synchrotron flux from the
MHD simulation output (namely, equipartition between the nonthermal particle
and magnetic energy densities, see Sect. 2), the synchrotron flux depends only
on magnetic field. However, there are two ways in which the magnetic field
structure can affect emission. On one hand, an overall increase in the mean
magnetic field strength will straightforwardly make the collision brighter. On
the other hand, synchrotron emission depends non-linearly on the magnetic
field strength (see Eq. 6). Thus, even if the mean magnetic field strength is
identical between 1D and 2D simulations, the presence of 2D structures in the
latter can still alter the resulting synchrotron flux. To investigate which of
these effects dominates in our case, we perform additional experiments. In a
selected 2D run with well-visible synchrotron emission enhancement (pressure-
corrugated case with $A_{P}=0.5$), we replace the magnetic field data with
various modifications that deactivate some of the mechanisms above. The
results are shown in Fig. 7. The blue curve therein shows the original result,
i.e., the 2D simulation with the original 2D magnetic field structure (c.f.,
Fig. 5). The orange curve corresponds to the same simulation where we replaced
its original 2D magnetic field with the 1D result (i.e., the 1D magnetic field
profile replicated along the $y$ direction to fill the 2D simulation box). As
this nullifies both of the mechanisms discussed – there is no change in mean
magnetic field and no 2D structure – this results in near-zero relative
difference333The remaining difference in flux is due to the 2D velocity
structure still present after replacing magnetic field data. The fluid
velocity affects transformation of magnetic fields from the observer frame
(where the magnetic fields are defined) to the fluid frame (where we calculate
synchrotron emission).. To obtain the green curve, we use the 2D magnetic
field result, but scale it so that the mean value matches that of the 1D
result. Thus, any enhancement in the mean value is removed, leaving only the
influence of dimensionality (2D structure). We see that the resulting curve
lies between the orange (no $F_{\rm syn}$ enhancement) and blue (the actual 2D
result). Thus, the 2D nature of the magnetic field spatial distribution
significantly contributes to the long-term synchrotron emission enhancement in
corrugated simulations, but magnetic field amplification (stronger in
corrugated than in 1D shocks), has comparable, if not slightly dominant,
influence.

###  3.2 Downstream magnetic field structure and polarization properties

Figure 8: Decay of the magnetic field component perpendicular to the shock
normal, $B_{\perp}$. Left: $B_{\perp}$ as a function of the distance from the
shock front, $x-x_{\rm sh}$. Solid lines show the vertically-averaged values
from the 2D simulations, dashed – from 1D runs. Right: slopes of the two power
laws fit to the measurements ($p_{1}$ – close to the shock front, $p_{2}$ –
further away), as a function of the corrugation aspect ratio. Open and full
symbols refer to the fits to 1D and 2D results, respectively. The top row
corresponds to density-corrugation simulations, while bottom row plots show
pressure-corrugation results. All curves and markers are color-coded by the
simulation amplitude, as indicated in the legend.

Corrugation of the shock front affects the structure of magnetic field in the
downstream medium (see Fig. 4 and Demidem et al. 2018). Information about the
geometry of magnetic field lines has been successfully extracted from
observations using polarimetry. In this section, we investigate how the
corrugated shell collisions affect the magnetic structure of the downstream
and the multiwavelength polarimetric properties. Synchrotron polarization is
orthogonal to the direction of magnetic field and is expected to be either
parallel or orthogonal to the jet direction for the cases of purely toroidal
and poloidal configurations, respectively. In this respect, the ratio of the
magnetic field components behind the shock front is of high importance.
Previously, the dominance of the parallel field component was assumed at the
shock front (Tavecchio et al., 2018), motivated by earlier studies of particle
energization possibilities at the shock fronts (Sironi & Spitkovsky, 2009;
Sironi et al., 2013). Tavecchio et al. (2018) find that the observed
polarization fractions of optical and X-ray emission in BL Lac jets can be
deduced from the slope of the power law fitting the decay of the shock-
perpendicular magnetic field $B_{\perp}$ (the component parallel to the shock
front, i.e., perpendicular to the shock normal) behind the shock.

In their analysis, Tavecchio et al. (2018) assume the pre-shock magnetic field
to be composed of a shock-parallel component (along the flow direction / the
shock normal) following a well-defined profile, and a randomly oriented shock-
perpendicular component. In this work, we consider a homogeneous shock-
perpendicular configuration. However, the post-shock magnetic field strength
is always locally dominated by the compressed shock-perpendicular component
(randomly oriented in the case of Tavecchio et al. 2018, see their fig. 3).
Thus, the post-shock region of our model can be viewed as a small region of
the post-shock domain of Tavecchio et al. (2018), where the magnetic field can
be treated as homogeneous.

Given these considerations, we can measure the decay slope of $B_{\perp}$ in
our simulations and discuss what it implies for the measured polarization in
light of the study by Tavecchio et al. (2018). After finding shock positions
using the local maxima of $|\partial B_{\perp}/\partial x|$, we plot the
$y$-averaged $B_{\perp}$ as a function of distance from the shock in Fig. 8.
Since each simulation contains two shocks (in the left and the right shell)
the average between the two is used.

We find that $B_{\perp}$ remains constant close to the shock and then drops,
approximately as a power law. This structure indicates that there is a
rarefaction trailing each shock, with a near-homogeneous fluid between them.
Due to the fact that the ambient medium between the shells (through which the
rarefaction travels) is adjusted between runs for density corrugation (see
Sect. 2.4), we observe much greater diversity of curve shapes in that case
(top row of Fig. 8), while runs with pressure corrugation (bottom row of Fig.
8) all behave similarly. In each case, the behavior of the matching 1D and 2D
runs is very similar, showing that the quantities linear in $B$, such as the
mean $B_{\perp}$ considered here, are well approximated by the solution to the
1D Riemann problem.

A (hard-)broken power law fit to $B_{\perp}(x-x_{\rm sh})$ yields the results
presented in the right column of Fig. 8. We see that the post-shock power law
slope remains within the range $[-0.50,-0.15]$ for density-corrugated
simulations and is well constrained at $\sim-0.14$ for the pressure-corrugated
ones. Tavecchio et al. (2018) find that the slope (in their notation, $-m$,
see their fig. 4) of $\sim-0.4$ results in high polarization fractions of both
the optical and X-ray band, with the latter somewhat higher than the former,
$\gtrsim 45\%$ vs $\sim 33\%$. We caution that these specific values would
only apply to the case of a single shell collision. Once a lot of events
overlap (with some collisions likely off-center, at different inclinations to
the jet axis), the expected total polarization fraction of jet emission should
be lower, and could be comparable to X-ray polarization recently detected by
IXPE (Liodakis et al., 2022). Nonetheless, unless another source of polarized
radiation contributes to the observed bands, the X-ray emission should remain
more strongly polarized than the optical band.

##  4 Discussion and conclusions

In this paper, we revise the viability of internal shock collisions to power
the emission observed from relativistic jets. Given sufficient presence of
high-energy particles, synchrotron emission in these environments is largely
dictated by the strength and orientation of magnetic fields (see discussion in
Sect. 2). At the base of the jet, these may be randomly, or even poloidally,
oriented (Vlahakis & Königl, 2003; Komissarov et al., 2009; Zdziarski et al.,
2022a). However, the magnetic fields parallel ($B_{||}$) and perpendicular
($B_{\perp}$) to the jet axis evolve differently with the distance from the
central engine (Blandford & Königl, 1979):

$\displaystyle B_{||}$ $\displaystyle\propto z^{-2}\textrm{ (magnetic flux
conservation)},$ (20) $\displaystyle B_{\perp}$ $\displaystyle\propto
z^{-1}\textrm{ (magnetic energy conservation)}.$

This causes the magnetic field to be predominantly toroidal in the radio-
emitting parts of the jet (as evidenced by observations of radio and optical
polarization of blazars, e.g., Gabuzda & Cawthorne 2000; Gabuzda et al. 2006;
Walker et al. 2018; Issaoun et al. 2022; Zhao et al. 2022).

Internal shocks are driven along the jet axis into this changing configuration
of magnetic fields. Close to the central engine, where $\vec{B}$ is randomly
or poloidally oriented, this may result in quasi-parallel shocks (with
$\vec{B}$ nearly parallel to the shock normal), which are known to be
efficient particle accelerators in the standard Diffusive Shock Acceleration
(DSA) mechanism (Sironi & Spitkovsky, 2009; Sironi et al., 2015a). Further
away, the internal shocks likely become nearly exclusively perpendicular (with
$\vec{B}$ perpendicular to the shock normal), the configuration that has long
been considered to be inefficient for particle acceleration (Begelman & Kirk,
1990; Sironi et al., 2013). This prompted efforts in considering alternative
particle energization mechanisms, such as magnetic reconnection (Sironi et
al., 2015b). It was recently found that perpendicular shocks may lead to
particle acceleration when their surface is corrugated (Demidem et al., 2022).
Such a situation is not unlikely. Turbulence and density inhomogeneities
within the jet can easily affect shock structure either directly or through
secondary instabilities (as has been shown for gamma-ray burst models, e.g.,
López-Cámara et al., 2016; Gottlieb et al., 2020).

While corrugation offers a way to produce suprathermal particles in quasi-
perpendicular magnetized shocks – a natural configuration for internal shocks
in jets – its net effect on the downstream structure and radiative properties
has not been investigated previously. Here, we assess whether the current
phenomenological approach to compute spectral and timing properties of the
internal shock model, which assumes one-dimentional collisions in their
emission calculations (e.g., Beloborodov, 2000; Malzac, 2014; Marino et al.,
2020), is still valid in the less trivial and arguably more realistic case of
corrugated shocks.

Assuming that particles attain a power-law energy distribution in both cases,
we find corrugated shocks to be up to $\sim 10\%$ brighter than their flat
counterparts. We acknowledge this to be a small correction in comparison to
effects of other approximations made in internal shock models. Thus, we
conclude that the existing internal shock models are entirely consistent with
the behavior of 2D corrugated shocks, at least in the simplified MHD approach
presented here. We caution that a detailed treatment of particle acceleration
by corrugated shocks (either through test particles or fully PIC simulations)
may well change this conclusion.

Multiwavelength polarization studies are well known to reveal the structure of
magnetic fields in jets (e.g. Gabuzda & Cawthorne, 2000; Gabuzda et al., 2006;
Issaoun et al., 2022). With the recent launch of the Imaging X-ray Polarimetry
Explorer (IXPE, Weisskopf et al. 2022), studies have for the first time been
able to reach the regions immediately behind the shock fronts. Our simulations
allow us to make predictions regarding the level of polarization of jet
synchrotron emission in the context of internal shock models.

From the point of view of a single shell collision, our results allow us to
apply the findings of Tavecchio et al. (2018), leading to the prediction that
a single two-shell collision in the region of toroidal magnetic field (far in
the jet) should result in low optical and high X-ray polarization fraction. At
the same time, a similar collision close to the central engine would occur at
a much higher value of $B_{||}$ (as discussed above), resulting in highly
polarized radiation in both spectral regions. However, the frequency of the
peak of the synchrotron emission flux depends on the distance from the central
engine ($\nu\mathchoice{\mathrel{\vbox{
\offinterlineskip\halign{\hfil$#$\cr\displaystyle\propto\cr\kern
2.0pt\cr\displaystyle\sim\cr\kern-2.0pt\cr}}}}{\mathrel{\vbox{
\offinterlineskip\halign{\hfil$#$\cr\textstyle\propto\cr\kern
2.0pt\cr\textstyle\sim\cr\kern-2.0pt\cr}}}}{\mathrel{\vbox{
\offinterlineskip\halign{\hfil$#$\cr\scriptstyle\propto\cr\kern
2.0pt\cr\scriptstyle\sim\cr\kern-2.0pt\cr}}}}{\mathrel{\vbox{
\offinterlineskip\halign{\hfil$#$\cr\scriptscriptstyle\propto\cr\kern
2.0pt\cr\scriptscriptstyle\sim\cr\kern-2.0pt\cr}}}}z^{-1}$, known as core-
shift, Lobanov 1998; Hirotani 2005), so regions close to the accretion disk
will mainly radiate in X-ray, while those further out will be seen in
optical/IR/radio bands. Thus, we can overall predict that the jet-integrated
synchrotron emission from an internal-shock-powered jet should exhibit a
highly-polarized X-ray emission and a weakly-polarized optical/IR emission. We
are excited to see these expectations systematically investigated with
multiwavelength polarimetric data.

Our work expands the 1D picture of shell collisions, typically used with
internal shock models, into 2D. A natural next step would be to investigate a
fully 3D system. Instabilities may behave differently in 2D than in 3D.
However, the only instability we observe in our runs is the Kelvin-Helmholtz
instability, which, by its nature, remains locally two-dimensional. As
reported by, e.g., Markwick et al. (2021) colliding flows may also exhibit
radiative instabilities which do evolve differently in 3D. Still, due to the
small radiative power of microquasar / AGN jets relative to their kinetic
power, these instabilities should not have significant dynamical consequences
on our systems of interest. Thus, we expect the downstream to behave similarly
in 3D as it does in 2D. The main finding of our 2D work is that magnetic field
is amplified in structured (as opposed to 1D) systems, leading to brighter
synchrotron emission. It is difficult to say whether this amplification would
be stronger or weaker in 3D configurations. However, a dynamo may operate in
3D in addition to the process we observe here, leading to further
amplification of magnetic fields. Thus, we expect the increase in synchrotron
brightness to be comparable or even larger in 3D than it is in 2D. Further
work in 3D would be needed to address whether the conjectures we list above
are correct.

This research was performed with computing resources (Beskow, Dardel)
generously provided by PDC Center for High Performance Computing at the KTH
Royal Instutute of Technology, allocated through the Swedish National
Infrastructure for Computing (SNIC). Nordita is supported in part by
NordForsk. C.D. acknowledges support from NSF grant AST 1903335 and NASA grant
NNX17AK55G. A.V. acknowledges support from the Academy of Finland grants
309308 and 322779.

## References

  * Bassi et al. (2020) Bassi, T., Malzac, J., Del Santo, M., et al. 2020, MNRAS, 494, 571, doi: 10.1093/mnras/staa739
  * Beckwith & Stone (2011) Beckwith, K., & Stone, J. M. 2011, ApJS, 193, 6, doi: 10.1088/0067-0049/193/1/6
  * Begelman & Kirk (1990) Begelman, M. C., & Kirk, J. G. 1990, ApJ, 353, 66, doi: 10.1086/168590
  * Beloborodov (2000) Beloborodov, A. M. 2000, ApJ, 539, L25, doi: 10.1086/312830
  * Blandford & Königl (1979) Blandford, R. D., & Königl, A. 1979, ApJ, 232, 34, doi: 10.1086/157262
  * Böttcher & Dermer (2010) Böttcher, M., & Dermer, C. D. 2010, ApJ, 711, 445, doi: 10.1088/0004-637X/711/1/445
  * Bright et al. (2020) Bright, J. S., Fender, R. P., Motta, S. E., et al. 2020, Nature Astronomy, 4, 697, doi: 10.1038/s41550-020-1023-5
  * Cerutti & Giacinti (2020) Cerutti, B., & Giacinti, G. 2020, A&A, 642, A123, doi: 10.1051/0004-6361/202038883
  * Connors et al. (2019) Connors, R. M. T., van Eijnatten, D., Markoff, S., et al. 2019, MNRAS, 485, 3696, doi: 10.1093/mnras/stz604
  * Corbel et al. (2013) Corbel, S., Coriat, M., Brocksopp, C., et al. 2013, MNRAS, 428, 2500, doi: 10.1093/mnras/sts215
  * Daigne & Mochkovitch (1998) Daigne, F., & Mochkovitch, R. 1998, MNRAS, 296, 275, doi: 10.1046/j.1365-8711.1998.01305.x
  * Davis & Tchekhovskoy (2020) Davis, S. W., & Tchekhovskoy, A. 2020, ARA&A, 58, 407, doi: 10.1146/annurev-astro-081817-051905
  * Demidem et al. (2018) Demidem, C., Lemoine, M., & Casse, F. 2018, MNRAS, 475, 2713, doi: 10.1093/mnras/stx3367
  * Demidem et al. (2022) Demidem, C., Nättilä, J., & Veledina, A. 2022, arXiv e-prints, arXiv:2212.06053.  https://arxiv.org/abs/2212.06053
  * Deng et al. (2015) Deng, W., Li, H., Zhang, B., & Li, S. 2015, ApJ, 805, 163, doi: 10.1088/0004-637X/805/2/163
  * Derishev et al. (2003) Derishev, E. V., Aharonian, F. A., Kocharovsky, V. V., & Kocharovsky, V. V. 2003, Phys. Rev. D, 68, 043003, doi: 10.1103/PhysRevD.68.043003
  * Drappeau et al. (2015) Drappeau, S., Malzac, J., Belmont, R., Gandhi, P., & Corbel, S. 2015, MNRAS, 447, 3832, doi: 10.1093/mnras/stu2711
  * Fan et al. (2004) Fan, Y. Z., Wei, D. M., & Zhang, B. 2004, MNRAS, 354, 1031, doi: 10.1111/j.1365-2966.2004.08263.x
  * Gabuzda & Cawthorne (2000) Gabuzda, D. C., & Cawthorne, T. V. 2000, MNRAS, 319, 1056, doi: 10.1046/j.1365-8711.2000.03881.x
  * Gabuzda et al. (2006) Gabuzda, D. C., Rastorgueva, E. A., Smith, P. S., & O’Sullivan, S. P. 2006, MNRAS, 369, 1596, doi: 10.1111/j.1365-2966.2006.10433.x
  * Gallo et al. (2005) Gallo, E., Fender, R., Kaiser, C., et al. 2005, Nature, 436, 819, doi: 10.1038/nature03879
  * Goodson et al. (1997) Goodson, A. P., Winglee, R. M., & Böhm, K.-H. 1997, ApJ, 489, 199, doi: 10.1086/304774
  * Gottlieb et al. (2020) Gottlieb, O., Levinson, A., & Nakar, E. 2020, MNRAS, 495, 570, doi: 10.1093/mnras/staa1216
  * Graff et al. (2008) Graff, P. B., Georganopoulos, M., Perlman, E. S., & Kazanas, D. 2008, ApJ, 689, 68, doi: 10.1086/592427
  * Granot (2012) Granot, J. 2012, MNRAS, 421, 2467, doi: 10.1111/j.1365-2966.2012.20474.x
  * Heida et al. (2017) Heida, M., Jonker, P. G., Torres, M. A. P., & Chiavassa, A. 2017, ApJ, 846, 132, doi: 10.3847/1538-4357/aa85df
  * Hirotani (2005) Hirotani, K. 2005, ApJ, 619, 73, doi: 10.1086/426497
  * Hynes et al. (2003) Hynes, R. I., Steeghs, D., Casares, J., Charles, P. A., & O’Brien, K. 2003, ApJ, 583, L95, doi: 10.1086/368108
  * Issaoun et al. (2022) Issaoun, S., Wielgus, M., Jorstad, S., et al. 2022, ApJ, 934, 145, doi: 10.3847/1538-4357/ac7a40
  * Jamil & Böttcher (2012) Jamil, O., & Böttcher, M. 2012, ApJ, 759, 45, doi: 10.1088/0004-637X/759/1/45
  * Jamil et al. (2010) Jamil, O., Fender, R. P., & Kaiser, C. R. 2010, MNRAS, 401, 394, doi: 10.1111/j.1365-2966.2009.15652.x
  * Joshi & Böttcher (2011) Joshi, M., & Böttcher, M. 2011, ApJ, 727, 21, doi: 10.1088/0004-637X/727/1/21
  * Kaiser et al. (2000) Kaiser, C. R., Sunyaev, R., & Spruit, H. C. 2000, A&A, 356, 975.  https://arxiv.org/abs/astro-ph/0001501
  * Kino et al. (2004) Kino, M., Mizuta, A., & Yamada, S. 2004, ApJ, 611, 1021, doi: 10.1086/422305
  * Kobayashi & Sari (2001) Kobayashi, S., & Sari, R. 2001, ApJ, 551, 934, doi: 10.1086/320249
  * Komissarov (2011) Komissarov, S. S. 2011, Mem. Soc. Astron. Italiana, 82, 95.  https://arxiv.org/abs/1006.2242
  * Komissarov et al. (2009) Komissarov, S. S., Vlahakis, N., Königl, A., & Barkov, M. V. 2009, MNRAS, 394, 1182, doi: 10.1111/j.1365-2966.2009.14410.x
  * Levine & Corbet (2006) Levine, A. M., & Corbet, R. 2006, The Astronomer’s Telegram, 940, 1
  * Ligorini et al. (2021) Ligorini, A., Niemiec, J., Kobzar, O., et al. 2021, MNRAS, 501, 4837, doi: 10.1093/mnras/staa3901
  * Liodakis et al. (2022) Liodakis, I., Marscher, A. P., Agudo, I., et al. 2022, arXiv e-prints, arXiv:2209.06227.  https://arxiv.org/abs/2209.06227
  * Lobanov (1998) Lobanov, A. P. 1998, A&A, 330, 79.  https://arxiv.org/abs/astro-ph/9712132
  * López-Cámara et al. (2016) López-Cámara, D., Lazzati, D., & Morsony, B. J. 2016, ApJ, 826, 180, doi: 10.3847/0004-637X/826/2/180
  * Lyutikov (2003) Lyutikov, M. 2003, New A Rev., 47, 513, doi: 10.1016/S1387-6473(03)00083-6
  * Lyutikov et al. (2012) Lyutikov, M., Balsara, D., & Matthews, C. 2012, MNRAS, 422, 3118, doi: 10.1111/j.1365-2966.2012.20831.x
  * MacFadyen & Woosley (1999) MacFadyen, A. I., & Woosley, S. E. 1999, ApJ, 524, 262, doi: 10.1086/307790
  * Maitra et al. (2009) Maitra, D., Markoff, S., Brocksopp, C., et al. 2009, MNRAS, 398, 1638, doi: 10.1111/j.1365-2966.2009.14896.x
  * Malzac (2013) Malzac, J. 2013, MNRAS, 429, L20, doi: 10.1093/mnrasl/sls017
  * Malzac (2014) —. 2014, MNRAS, 443, 299, doi: 10.1093/mnras/stu1144
  * Malzac et al. (2018) Malzac, J., Kalamkar, M., Vincentelli, F., et al. 2018, MNRAS, 480, 2054, doi: 10.1093/mnras/sty2006
  * Marino et al. (2020) Marino, A., Malzac, J., Del Santo, M., et al. 2020, MNRAS, 498, 3351, doi: 10.1093/mnras/staa2570
  * Markoff et al. (2005) Markoff, S., Nowak, M. A., & Wilms, J. 2005, ApJ, 635, 1203, doi: 10.1086/497628
  * Markwick et al. (2021) Markwick, R. N., Frank, A., Carroll-Nellenback, J., et al. 2021, MNRAS, 508, 2266, doi: 10.1093/mnras/stab2577
  * Matthews et al. (2020) Matthews, J. H., Bell, A. R., & Blundell, K. M. 2020, New A Rev., 89, 101543, doi: 10.1016/j.newar.2020.101543
  * Medvedev & Loeb (1999) Medvedev, M. V., & Loeb, A. 1999, ApJ, 526, 697, doi: 10.1086/308038
  * Meyer et al. (2015) Meyer, E. T., Georganopoulos, M., Sparks, W. B., et al. 2015, Nature, 521, 495, doi: 10.1038/nature14481
  * Mignone et al. (2009) Mignone, A., Ugliano, M., & Bodo, G. 2009, MNRAS, 393, 1141, doi: 10.1111/j.1365-2966.2008.14221.x
  * Mimica & Aloy (2010) Mimica, P., & Aloy, M. A. 2010, MNRAS, 401, 525, doi: 10.1111/j.1365-2966.2009.15669.x
  * Mimica & Aloy (2012) —. 2012, MNRAS, 421, 2635, doi: 10.1111/j.1365-2966.2012.20495.x
  * Mimica et al. (2007) Mimica, P., Aloy, M. A., & Müller, E. 2007, A&A, 466, 93, doi: 10.1051/0004-6361:20066811
  * Mimica et al. (2004) Mimica, P., Aloy, M. A., Müller, E., & Brinkmann, W. 2004, A&A, 418, 947, doi: 10.1051/0004-6361:20034261
  * Mimica et al. (2005) —. 2005, A&A, 441, 103, doi: 10.1051/0004-6361:20053218
  * Mirabel & Rodríguez (1999) Mirabel, I. F., & Rodríguez, L. F. 1999, ARA&A, 37, 409, doi: 10.1146/annurev.astro.37.1.409
  * Péault et al. (2019) Péault, M., Malzac, J., Coriat, M., et al. 2019, MNRAS, 482, 2447, doi: 10.1093/mnras/sty2796
  * Pe’er (2014) Pe’er, A. 2014, Space Sci. Rev., 183, 371, doi: 10.1007/s11214-013-0001-y
  * Pe’er et al. (2017) Pe’er, A., Long, K., & Casella, P. 2017, ApJ, 846, 54, doi: 10.3847/1538-4357/aa80df
  * Rees (1978) Rees, M. J. 1978, MNRAS, 184, 61P, doi: 10.1093/mnras/184.1.61P
  * Rees & Meszaros (1994) Rees, M. J., & Meszaros, P. 1994, ApJ, 430, L93, doi: 10.1086/187446
  * Rudolph et al. (2020) Rudolph, A., Heinze, J., Fedynitch, A., & Winter, W. 2020, ApJ, 893, 72, doi: 10.3847/1538-4357/ab7ea7
  * Rueda-Becerril et al. (2014) Rueda-Becerril, J. M., Mimica, P., & Aloy, M. A. 2014, MNRAS, 438, 1856, doi: 10.1093/mnras/stt2335
  * Rueda-Becerril et al. (2015) —. 2015, arXiv e-prints, arXiv:1502.07882.  https://arxiv.org/abs/1502.07882
  * Rybicki & Lightman (1979) Rybicki, G. B., & Lightman, A. P. 1979, Radiative processes in astrophysics
  * Shidatsu et al. (2011) Shidatsu, M., Ueda, Y., Tazaki, F., et al. 2011, PASJ, 63, S785, doi: 10.1093/pasj/63.sp3.S785
  * Sironi et al. (2015a) Sironi, L., Keshet, U., & Lemoine, M. 2015a, Space Sci. Rev., 191, 519, doi: 10.1007/s11214-015-0181-8
  * Sironi et al. (2015b) Sironi, L., Petropoulou, M., & Giannios, D. 2015b, MNRAS, 450, 183, doi: 10.1093/mnras/stv641
  * Sironi & Spitkovsky (2009) Sironi, L., & Spitkovsky, A. 2009, ApJ, 698, 1523, doi: 10.1088/0004-637X/698/2/1523
  * Sironi & Spitkovsky (2011) —. 2011, ApJ, 726, 75, doi: 10.1088/0004-637X/726/2/75
  * Sironi et al. (2013) Sironi, L., Spitkovsky, A., & Arons, J. 2013, ApJ, 771, 54, doi: 10.1088/0004-637X/771/1/54
  * Spada et al. (2001) Spada, M., Ghisellini, G., Lazzati, D., & Celotti, A. 2001, MNRAS, 325, 1559, doi: 10.1046/j.1365-8711.2001.04557.x
  * Stern & Poutanen (2008) Stern, B. E., & Poutanen, J. 2008, MNRAS, 383, 1695, doi: 10.1111/j.1365-2966.2007.12706.x
  * Stone & Gardiner (2009) Stone, J. M., & Gardiner, T. 2009, New A, 14, 139, doi: 10.1016/j.newast.2008.06.003
  * Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137, doi: 10.1086/588755
  * Tavecchio et al. (2018) Tavecchio, F., Landoni, M., Sironi, L., & Coppi, P. 2018, MNRAS, 480, 2872, doi: 10.1093/mnras/sty1491
  * Tchekhovskoy et al. (2009) Tchekhovskoy, A., McKinney, J. C., & Narayan, R. 2009, ApJ, 699, 1789, doi: 10.1088/0004-637X/699/2/1789
  * Tetarenko et al. (2019) Tetarenko, A. J., Casella, P., Miller-Jones, J. C. A., et al. 2019, MNRAS, 484, 2987, doi: 10.1093/mnras/stz165
  * Toro (1999) Toro, E. F. 1999, Riemann Solvers and Numerical Methods for Fluid Dynamics (2nd ed.; Berlin: Springer)
  * Vincentelli et al. (2019) Vincentelli, F. M., Casella, P., Petrucci, P., et al. 2019, ApJ, 887, L19, doi: 10.3847/2041-8213/ab5860
  * Vlahakis & Königl (2003) Vlahakis, N., & Königl, A. 2003, ApJ, 596, 1080, doi: 10.1086/378226
  * Walker et al. (2018) Walker, R. C., Hardee, P. E., Davies, F. B., Ly, C., & Junor, W. 2018, ApJ, 855, 128, doi: 10.3847/1538-4357/aaafcc
  * Weisskopf et al. (2022) Weisskopf, M. C., Soffitta, P., Baldini, L., et al. 2022, Journal of Astronomical Telescopes, Instruments, and Systems, 8, 026002, doi: 10.1117/1.JATIS.8.2.026002
  * Zdziarski et al. (2012) Zdziarski, A. A., Lubiński, P., & Sikora, M. 2012, MNRAS, 423, 663, doi: 10.1111/j.1365-2966.2012.20903.x
  * Zdziarski et al. (2015) Zdziarski, A. A., Sikora, M., Pjanka, P., & Tchekhovskoy, A. 2015, MNRAS, 451, 927, doi: 10.1093/mnras/stv986
  * Zdziarski et al. (2022a) Zdziarski, A. A., Stawarz, Ł., Sikora, M., & Nalewajko, K. 2022a, MNRAS, 515, L17, doi: 10.1093/mnrasl/slac060
  * Zdziarski et al. (2022b) Zdziarski, A. A., Tetarenko, A. J., & Sikora, M. 2022b, ApJ, 925, 189, doi: 10.3847/1538-4357/ac38a9
  * Zhao et al. (2022) Zhao, G.-Y., Gómez, J. L., Fuentes, A., et al. 2022, ApJ, 932, 72, doi: 10.3847/1538-4357/ac6b9c
  * Zhdankin et al. (2017) Zhdankin, V., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2017, Phys. Rev. Lett., 118, 055103, doi: 10.1103/PhysRevLett.118.055103

##  Appendix A Parameter range estimation for GX 339$-$4 environment

To obtain plasma parameters close to what is expected in GX 339$-$4, we use
existing theoretical models of the binary from the literature to provide lower
and upper constraints on the values of interest, and then select values within
the resulting ranges.

##### Lower limits:

Maitra et al. (2009) calculates the particle density of leptons at the jet
ejection zone of GX 339$-$4 to be:

$n_{\rm e}(r_{\rm ej})\sim 4\times[10^{13.4},10^{15}]\,\textrm{cm}^{-3}\gtrsim
10^{14}\,\textrm{cm}^{-3},$ (A1)

where the factor of $4$ comes from transforming to post-shock density. The
location of the ejection region is constrained within (their Table 2):

$z_{\rm ej}\sim 33-114r_{\rm g}.$ (A2)

Here, we will assume the jet to be roughly conical between the ejection zone
and the IR zone of Shidatsu et al. (2011) at $\sim 10^{4}r_{\rm g}$. As the
jet is likely convex while accelerating, this will likely underestimate the
local particle number density:

$n_{\rm e}(z)\propto z^{-2},$ (A3) $n_{\rm e}(z_{\rm IR})\gtrsim 1.3\times
10^{10}\,\textrm{cm}^{-3}.$ (A4)

This gives us the local mass density:

$\rho(z_{\rm IR})\gtrsim n_{\rm e}m_{\rm p}\sim 4.2\times 10^{-14}\textrm{
g~{}cm}^{-3}$ (A5)

and magnetization:

$\sigma_{\rm j}(z_{\rm IR})\simeq\frac{B^{2}/8\pi}{\rho(z_{\rm
IR})c^{2}}\lesssim 2.6.$ (A6)

We can also calculate the plasma skin depth:

$l_{\rm s}=\frac{c}{\omega_{pe}}=c\sqrt{\frac{m_{e}}{4\pi n_{\rm
e}e^{2}}}\lesssim 3.4\textrm{ cm}.$ (A7)

##### Upper limits:

As argued by Maitra et al. (2009), the total jet power $P_{\rm j}$ is larger
than bulk kinetic energy advection in protons by a factor of $\sim 10$, as it
also contains the energy later used to reheat the particle distribution (by,
e.g., internal shocks) (see discussion in their Sect. 2.1 regarding the “input
power parameter” $N_{\rm j}$ and its relation to total jet power of Markoff et
al. 2005). Interestingly, their $N_{\rm j}$ (see their Table 2) is actually
very close in value to $P_{\rm j}$ from Malzac et al. (2018) (however, as
these are physically different models, this observation need not be
necessarily meaningful). Given these findings, as well as results regarding
jet power from Gallo et al. (2005) and Bright et al. (2020), lowering the jet
power $P_{\rm j}$ of Malzac et al. (2018) by a factor of 10 can be used to
estimate the local plasma parameters. Thus, in this paragraph, we will take
$P_{\rm j}\sim 10^{37}\textrm{ erg}\textrm{s}^{-1}$.

Assuming that jet power is fully dominated by the kinetic energy of protons
(i.e., overestimating jet power):

$P_{\rm j}\lesssim(\Gamma_{\rm j}-1)\rho c^{2}\times v_{\rm
j}\times\pi(\theta_{\rm j}z)^{2},$ (A8) $\rho(z)\lesssim\frac{P_{\rm
j}}{(\Gamma_{\rm j}-1)v_{\rm j}c^{2}\times\pi(\theta_{\rm j}z)^{2}},$ (A9)
$\sigma_{\rm j}(z)\gtrsim\frac{B^{2}(\Gamma_{\rm j}-1)v_{\rm
j}\times\pi(\theta_{\rm j}z)^{2}}{8\pi P_{\rm j}}.$ (A10)

This, for the location of IR core from Shidatsu et al. (2011), gives us:

$\rho(z_{\rm IR})\lesssim 4.3\times 10^{-13}\textrm{ g~{}cm}^{-3},$ (A11)
$\sigma_{\rm j}(z_{\rm IR})\gtrsim 0.25.$ (A12)

Note that this is given in the frame of reference comoving with the jet at its
bulk Lorentz factor, $\Gamma_{\rm j}\sim 2$.

Assuming charge neutrality, we can obtain electron number density in this
frame from the mass density above:

$n_{\rm e}\simeq\rho(z)/m_{\rm p}\gtrsim 2.6\times 10^{11}\,\textrm{cm}^{-3}.$
(A13)

This gives a plasma skin depth of:

$l_{\rm s}=\frac{c}{\omega_{pe}}=c\sqrt{\frac{m_{e}}{4\pi n_{\rm
e}e^{2}}}\lesssim 1.0\textrm{ cm}.$ (A14)

##### Combined constraints

Based on the considerations above, we can constrain the fluid parameters
within GX 339$-$4 jet to be within the following ranges:

$\rho\in[4.2\times 10^{-14},4.3\times 10^{-13}]\textrm{ g~{}cm}^{-3},$ (A15)
$n_{\rm e}\in[1.3\times 10^{10},2.6\times 10^{11}]\textrm{ cm}^{-3},$ (A16)
$\beta\sim 1,$ (A17) $\sigma\in[2.6,0.25],$ (A18) $l_{\rm
s}\in[3.4,1.0]\textrm{ cm},$ (A19)

where $\rho$ denotes density, $n_{\rm e}$ – electron density, $\beta$ – plasma
beta (defined as the ratio of particle to magnetic field energy density),
$\sigma$ – magnetization, and $l_{\rm s}$ is the plasma skin depth. Note that
we have $\sim 10^{8}-10^{9}$ of local skin depths per jet radius at our region
of interest at $z_{\rm IR}\sim 10^{4}r_{\rm g}$.