# Exploring the role of composition and mass-loading on the properties of
hadronic jets

D. Kantzas,1,2 S. Markoff$,^{1,2}$ M. Lucchini,3 C. Ceccobello4 and K.
Chatterjee5
1Anton Pannekoek Institute for Astronomy (API), University of Amsterdam,
Science Park 904, 1098 XH Amsterdam, the Netherlands
2GRavitation AstroParticle Physics Amsterdam (GRAPPA), University of
Amsterdam, Science Park 904, 1098 XH Amsterdam, the Netherlands
3MIT Kavli Institute for Astrophysics and Space Research, Massachusetts
Institute of Technology, Cambridge, MA 02139, the USA
4Department of Space, Earth and Environment, Chalmers University of
Technology, Onsala Space Observatory, 439 92 Onsala, Sweden
5Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA
02138, the USA
E-mail<EMAIL_ADDRESS>at: LAPTh, CNRS, USMB, F-74940 Annecy,
France

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Astrophysical jets are relativistic outflows that remain collimated for
remarkably many orders of magnitude. Despite decades of research, the origin
of cosmic rays (CRs) remains unclear, but jets launched by both supermassive
black holes in the centre of galaxies and stellar-mass black holes harboured
in X-ray binaries (BHXBs) are among the candidate sources for CR acceleration.
When CRs accelerate in astrophysical jets, they initiate particle cascades
that form $\gamma$-rays and neutrinos. In the so-called hadronic scenario, the
population of accelerated CRs requires a significant amount of energy to
properly explain the spectral constraints similarly to a purely leptonic
scenario. The amount of energy required often exceeds the Eddington limit, or
even the total energy available within the jets. The exact energy source for
the accelerated protons is unclear, but due to energy conservation along the
jets, it is believed to come from the jet itself via transfer of energy from
the magnetic fields, or kinetic energy from the outflow. To address this
hadronic energy issue and to self-consistently evolve the energy flux along
the flows, we explore a novel treatment for including hadronic content, in
which instabilities along the jet/wind border play a critical role. We discuss
the impact of the different jet composition on the jet dynamics for a pair
dominated and an electron-proton jet, and consequently the emitted spectrum,
accounting for both leptonic and hadronic processes. Finally, we discuss the
implications of this mass-loading scenario to address the proton energy issue.

###### keywords:

acceleration of particles – stars: jets – galaxies: jets

††pubyear: 2022††pagerange: Exploring the role of composition and mass-loading
on the properties of hadronic jets–E

##  1 Introduction

Accreting black holes can efficiently launch relativistic outflows, known as
astrophysical jets, by converting gravitational energy to kinetic energy.
Large-scale jets launched by supermassive black holes (SMBH) share some common
physical laws to the small-scale jets launched by stellar-mass black holes in
X-ray binaries (BHXBs; Heinz & Sunyaev 2003; Merloni et al. 2003; Falcke et
al. 2004), and hence black hole jets appear to be scale invariant in some of
their properties. For example, SMBHs with masses of the order of $\sim
10^{6}-10^{9}\,M_{\rm\odot}$ power jets that remain collimated up to Mpc
scales (Waggett et al., 1977), whereas BHXBs with mass of the order of a few
solar masses display jets that remain collimated up to sub-pc scales (Mirabel
& Rodriguez, 1994). Galactic BHXBs are of particular importance because they
transition between different jetted and non-jetted states over human-like
timescales, giving us the chance to understand plasma evolution in extreme
conditions and better probe jet physics (see, e.g., Markoff et al., 2001;
Markoff et al., 2003; Markoff et al., 2005; Reig et al., 2003; Giannios et
al., 2004; Maitra et al., 2009; Vila & Romero, 2010; Zdziarski et al., 2014;
Connors et al., 2019; Lucchini et al., 2021).

The exact physical mechanism responsible for jet launching is not clear yet.
On one hand, the Blandford-Znajek mechanism (Blandford & Znajek, 1977)
describes a way to extract the rotational energy of a spinning black hole and
power relativistic jets that can be pair-plasma dominated (see, e.g.,
Broderick & Tchekhovskoy, 2015; Parfrey et al., 2019) On the other hand,
magnetic fields anchored in the accretion disc can launch baryon/proton/ion-
dominated jets via the Blandford-Payne mechanism (Blandford & Payne, 1982).
The difference in jet composition from the two launching mechanisms would have
an important impact on the interpretation of the spectral energy distribution
(SED) observed from such black hole systems as well as the consideration of
relativistic jets as candidate sources of cosmic rays (CRs).

CRs are charged particles that exhibit a large range of energies going up to
ultra-high energies of the order of $10^{20}\,\rm eV$ (The Pierre Auger
Observatory et al., 2017; Abbasi et al., 2020). The detected CR spectrum shows
two very prominent features, known as the “knee”and the “ankle” where the
spectrum steepens and hardens, respectively. The “knee” is observed at
$10^{15}\,\rm eV$ (PeV) and is likely to be the maximum energy that CR protons
accelerated in Galactic sources can reach, but the identification of these
particular sources remains a mystery despite the decades of studies. The
“ankle”, located at $\sim 10^{18}\,\rm eV$ (EeV), is where extragalactic
sources are thought to start dominating the spectrum. The exact CR composition
is not clear and strongly depends on the particle energy. GeV CRs primarily
comprise of protons ($\sim 99$ per cent; Shikaze et al. 2007), with electrons
and positrons mainly contributing to the rest of the spectrum. It is likely
that heavier elements/ions accelerated in Galactic sources start dominating
the CR spectrum between the “knee” and the “ankle” (Aloisio et al., 2012),
beyond which the composition is unclear (Abbasi et al., 2019; Yushkov et al.,
2019; Corstanje et al., 2021).

Similar to large-scale jets of active galactic nuclei (AGN), which are among
the dominant candidate sources of the extragalactic CRs (Protheroe & Kazanas,
1983), recent studies suggest the small-scale jets of BHXBs as potential CR
acceleration sites (Romero et al., 2003; Fender et al., 2005; Cooper et al.,
2020). There are currently only a few tens of Galactic BHXBs detected in the
Milky Way (Tetarenko et al., 2016), but population-synthesis simulations (see,
e.g., Olejak et al., 2020) suggest that a few thousand black holes likely
reside in the Galactic disc, in agreement with the recent X-ray observations
of the Galactic centre by Hailey et al. (2018) and Mori et al. (2021). Based
on such observations, Cooper et al. (2020) proposed that a few thousand BHXBs
are capable of contributing to the observed CR spectrum above the “knee”.

Whether or not BHXBs jets can indeed accelerate CRs up to the “knee”, and AGN
jets beyond the “ankle”, strongly depends on two further issues: (1) can
astrophysical jets, in general, accelerate particles to high energies, and (2)
are astrophysical jets actually comprised of protons and/or heavier elements?
On the former, observations of non-thermal emission from radio bands (see,
e.g., Lister et al., 2016) up to GeV/TeV $\gamma$-rays from both SMBHs (see,
e.g., Lister et al., 2009) and BHXBs (see, e.g., Zanin et al., 2016), suggest
that both classes of jets can efficiently accelerate particles. Numerous
numerical studies, moreover, suggest that jets can indeed be viable sites of
particle acceleration either via shocks (Hillas, 1984), or via magnetic
reconnection (Drenkhahn & Spruit, 2002; Guo et al., 2014; Sironi & Spitkovsky,
2014; Matthews et al., 2020).

The jet composition however remains an open question. The two different
proposed launching mechanisms mentioned above yield an entirely different jet
content at the base that significantly alters not only the jet dynamics, but
the emitted spectrum as well (Petropoulou et al., 2019). A pair-dominated jet
would allow only for leptonic processes, such as synchrotron and inverse
Compton scattering (ICS; Blumenthal & Gould 1970). A leptonic plus hadronic
jet, on the other hand, allows for further non-thermal processes, when
inelastic collisions occur between the accelerated protons and the cold flow
or radiation (e.g., Mannheim, 1993; Rachen & Biermann, 1993; Mannheim &
Schlickeiser, 1994; Rachen & Mészáros, 1998). Such hadronic processes can lead
to the production of astrophysical neutrinos, but usually require a much
larger jet energy budget than the leptonic ones, sometimes requiring super-
Eddington jet powers (Böttcher et al., 2013; Liodakis & Petropoulou, 2020).
Such super-Eddington powers challenge the accretion paradigm (Zdziarski &
Böttcher, 2015), but they still seem feasible for relativistic AGN jets
(Ghisellini et al., 2014).

Several BHXB jets, such as the peculiar case of SS433 or the prototypical
Cygnus X–1, show evidence of baryonic jet content (Fabrika 2004 and Gallo et
al. 2005; Heinz 2006, respectively). Both the compact objects of SS433 and
Cygnus X–1 are accompanied by a high mass donor star that may be the source of
the heavy composition through its stellar wind. There is evidence of baryon-
loaded jets though, even in the case of a low-mass companion, such as the
black hole candidate 4U 1630–47, based on iron emission lines (Díaz Trigo et
al., 2013). The cases of MAXI J1820+070 (Tetarenko et al., 2021; Zdziarski et
al., 2022a), MAXI J1836-194 (Lucchini et al., 2021), XTE J1752–223, MAXI
J1659–152, and XTE J1650–500 (Cao et al., 2021) on the other hand, favour a
jet composition of the order of a few to a few tens of pairs per proton based
on energetic arguments.

The composition is also difficult to constrain in extragalactic jets. Circular
polarisation measurements indicate that the jets of the blazar 3C 279 are
pair-dominated (Liodakis et al., 2021), and energetic arguments of the radio
galaxy 3C 120 are consistent with a pair-dominated jet (Zdziarski et al.,
2022b). Celotti & Fabian (1993), on the other hand, based on very-large
baseline interferometry and spectral arguments for numerous sources, support
an electron-proton plasma. The blazar TXS 0506+056, finally, due to the
correlation with the high-energy neutrino IceCube-170922A, supports a baryon
content in its jets as well (Aartsen et al., 2018).

Currently, the state-of-the-art to model jet launching and dynamics in a more
a priori way are high-resolution simulations that solve the magneto
hydrodynamic equations in the general relativistic regime (GRMHD). Such
simulations have furthered our understanding of the accretion-launching
paradigm and have shown that a Poynting flux dominated outflow can convert a
significant amount of its initial magnetic energy into kinetic energy to
accelerate the bulk flow (McKinney, 2006; Komissarov et al., 2007;
Tchekhovskoy et al., 2008, 2009; Komissarov et al., 2009). The same
simulations, have established that the accretion disc can significantly impact
the spatial evolution of the jets not only at $r_{g}$-scale distances
($r_{g}=GM_{\rm bh}/\rm c^{2}$, where $M_{\rm bh}$ is the mass of the black
hole), but also further out. In particular, Chatterjee et al. (2019, hereafter
CLTM19) performed a series of high-resolution GRMHD simulations of strongly
magnetised systems to better understand the loading of jets with matter from
the wind of the accretion disc. When the jets propagate in a medium, pinch
instabilities can occur in the interface between the jet and the ambient
medium to give rise to eddies that eventually allow for matter to entrain the
jet (Eichler 1993; Spruit et al. 1997; Begelman 1998; Giannios & Spruit 2006;
CLTM19; Sironi et al. 2021). Such mass entrainment can significantly affect
the jet kinematics and hence the non-thermal emission.

Such GRMHD simulations, though, usually make the ideal gas assumption and
therefore cannot capture dissipative processes like particle acceleration
self-consistently. Kinetic simulations of particles-in-cell (PIC), on the
other hand, calculate the trajectories of individual particles based on first
principles, allowing for a more detailed and comprehensive understanding of
the relativistic outflows. Both GRMHD and PIC simulations, however, are very
computational expensive, and they cannot easily be compared to observations
through statistical methods that explore the full parameter phase space.

In this work, we develop a new treatment for incorporating mass-loading and
thus evolving compositions in jets, and apply it to a multi-zone jet model.
This treatment is inspired by recent GRMHD simulations such as CLTM19, to
explore jet composition and its impact on the total jet power as well as its
electromagnetic emission. In particular, we build on the multi-zone jet model
developed by Markoff et al. (2005) that relies on the pioneering ideas of
Blandford & Königl (1979), Hjellming & Johnston (1988), and Falcke & Biermann
(1995). After many developments, the latest version of the model is BHJet
(Lucchini et al., 2022), a multi-zone jet model that better connects the jet
acceleration and jet physical quantities to the radiative output. For the
first time, we connect the physically motivated model BHJet with hadronic
acceleration, accounting for self-consistent energy conservation. We further
present HadJet, a multi-zone, lepto-hadronic, mass-loaded jet model. In this
work, we discuss the main physical properties of both models and how HadJet
can be used to address the jet-power crisis of lepto-hadronic models.

The paper is structured as follows. In Section 2 we describe the semi-
analytical calculations for the magnetically accelerated jet accounting for
both leptonic and hadronic acceleration and radiative processes. We present
the results of the above jet model in Section 3. In Section 4, we describe the
details of the mass-loaded jet model (HadJet) and present the results in
Section 5. Finally, in Section 6 we discuss the implication of our new models
on the proton power issue and conclude in Section 7.

##  2 Magnetically accelerated steady-state jets

We assume two initially cold, Poynting flux dominated jets of either leptonic
or lepto-hadronic content, that accelerate up to some maximum velocity because
of magnetic energy dissipation (Vlahakis & Konigl, 2003; McKinney, 2006;
Komissarov et al., 2007). At the region where the bulk velocity reaches the
maximum value (acceleration region henceforth, denoted by $z_{\rm acc}$), we
further assume that energy is also dissipated to accelerate particles to non-
thermal energies (Blandford & Rees, 1974; Begelman et al., 1984). With our
formalism, we cannot capture whether the magnetic energy dissipates
immediately to particle acceleration (as in the case of magnetic reconnection)
or if magnetic energy dissipates to kinetic energy first and this extra
kinetic energy dissipates to particle acceleration through shocks (Bogovalov &
Tsinganos, 2005). We assume instead that the total energy of the jet is
conserved at the particle acceleration region. From this point outwards along
the jets, we assume a constant particle acceleration rate and discuss below
how this assumption affects the evolution of both the jet velocity and
magnetic field. In Table 1, we define all the parameters and their fiducial
values (if applicable) that we use in this section.

###  2.1 Jet dynamical properties

Based on both semi-analytical and numerical calculations, the bulk jet Lorentz
factor $\gamma$ is expected to scale approximately as $z^{1/2}$, where $z$ is
the distance along the jet (Beskin & Nokhrina, 2006; McKinney, 2006). We
parametrise the jet Lorentz factor as Lucchini et al. (2018) (and see also
Potter & Cotter 2012)

$\gamma(z\leq z_{\rm acc})=\gamma_{0}+(\gamma_{\rm
acc}-\gamma_{0})\dfrac{z^{1/2}-z_{0}^{1/2}}{z_{\rm acc}^{1/2}-z_{0}^{1/2}},$
(1)

where $\gamma_{0}$ is the initial Lorentz factor at the jet base and $z_{0}$
is the distance of the jet base from the black hole and $\gamma$ acc is the
maximum bulk Lorentz factor at $z_{\rm diss}$. We assume that the jets launch
initially with the speed of sound, which for a relativistic flow with
adiabatic index 4/3 is equal to 0.43 c, or $\gamma_{0}=1.11$ (Crumley et al.,
2017).

The jets are thus set to be initially parabolic while they accelerate and
become conical when they achieve $\gamma_{\rm acc}$ (Komissarov et al., 2009).
We express the cross-sectional radius of the jet along the jet axis as

$r=r_{0}+(z-z_{0})\tan(\theta),$ (2)

where $r_{0}$ is the radius of the jet base and $\theta$ is the opening angle
of the jets. Based on very long baseline interferometry observations and the
Monitoring of jets in AGN with VLBA Experiments (MOJAVE; Pushkarev et al.
2009, 2017, see, e.g.,, we set the jet opening angle to be

$\theta=\dfrac{0.15}{\gamma}.$ (3)

While the number of particles along the jet is conserved, we express the
number density of leptons as

$n=n_{0}\left(\dfrac{\gamma\beta}{\gamma_{0}\beta_{0}}\right)^{-1}\left(\dfrac{r}{r_{0}}\right)^{-2},$
(4)

where $\beta$ is the jet velocity normalized to the speed of light and $n_{0}$
is the initial number density. We calculate $n_{0}$ by the power $L_{\rm jet}$
injected at the jet base in the comoving frame

$L_{\rm jet}=2\beta_{0}\gamma_{0}{\rm c\pi}r_{0}^{2}\omega_{0}$ (5)

where we account for two identical jets (hence the factor of 2), and $n_{0}$
depends on $L_{\rm jet}$ and the initial conditions of the jet base as written
out below. We write the jet enthalpy $\omega$ as (Falcke & Biermann, 1995;
Crumley et al., 2017)

$\omega=\rho{\rm c^{2}}+U_{j}+P_{j}=\rho c^{2}+U_{\rm p}+P_{\rm p}+U_{\rm
e}+P_{\rm e}+U_{B}+P_{B},$ (6)

where $U_{j}=U_{\rm p}+U_{\rm e}+U_{B}$ is the total internal jet energy
density and $P_{j}=P_{\rm p}+P_{\rm e}+P_{B}$ is the total jet pressure. In
the above equation, $\rho$ is the jet mass density

$\rho=n_{\rm p}\,{\rm m_{p}}+n_{\rm e}\,\rm m_{e}.$ (7)

We express the number of protons in terms of the number of leptons as $n_{\rm
p}=n_{\rm e}/\eta_{e}$, where $n_{\rm e/p}$ is the number density of
leptons/protons, respectively, and $\eta_{e}\geq 1$ is a free parameter that
remains constant unless the jets are mass-loaded (see below).

For an ideal gas, we can write the pressure terms as

$P_{\rm e,p}=(\Gamma_{\rm e,p}-1)\,U_{\rm e,p},$ (8)

where $\Gamma_{\rm e,p}$ is the adiabatic index. For the rest of the paper, we
assume a relativistic pair content ($\Gamma_{\rm e}=4/3$) at the jet base and
a cold proton population ($\Gamma_{\rm p}=5/3$) until the particle
acceleration region (see below). For the pair temperatures we are interested
in this work, the flow remains cold even if is dominated by pairs at the base.
For $U_{B}=P_{B}=B^{2}/8\pi$, we write the jet enthalpy as

$\omega=\rho{\rm c^{2}}+\Gamma_{\rm p}U_{\rm p}+\Gamma_{\rm e}U_{\rm
e}+\dfrac{B^{2}}{4\pi}.$ (9)

We define the specific enthalpy of the gas as

$h=\dfrac{U_{g}+P_{g}}{\rho{\rm{c^{2}}}}=\dfrac{\Gamma_{\rm p}U_{\rm
p}+\Gamma_{\rm e}U_{\rm e}}{\rho{\rm{c^{2}}}}$ (10)

where we used equation (8). We calculate $U_{\rm e,p}$ by computing the
integral

$U_{\rm e,p}=\int\dfrac{dn_{\rm e,p}}{d\varepsilon_{\rm e,p}}\varepsilon_{\rm
e,p}\rm m_{e,p}c^{2}d\varepsilon_{\rm e,p}.$ (11)

where $\varepsilon_{\rm e,p}$ is the Lorentz factor of the particles, but we
can also express the internal energy density in terms of the average total
energy of the particles

$U_{\rm e,p}\simeq(\langle\varepsilon_{\rm e,p}\rangle-1)\,n_{\rm e,p}\rm
m_{e,p}c^{2},$ (12)

where $\langle\varepsilon_{\rm e,p}\rangle$ is the average Lorentz factor of
the pairs/protons of the jet segment (see below for calculation). This
equation is more convenient than equation (11) for the following discussion,
however we note that it might not be accurate enough if a significant fraction
of the leptons accelerate to non-thermal energies, in particular in a hard
power law with slope $<2$.

A useful parameter to characterise the jets is the magnetisation. We define
the magnetisation of a flow as the Poynting flux over the total energy flux
(Nokhrina et al., 2015)

$\begin{split}&\sigma=\dfrac{B^{2}}{4\pi\left(\rho\rm{c^{2}}+U_{g}+P_{g}\right)}\Rightarrow\\\
&\sigma=\dfrac{B^{2}}{4\pi\rho{\rm c^{2}}\left(1+h\right)}.\end{split}$ (13)

When the flow is cold ($h\ll 1$), the above definition reduces to the well-
known expression of

$\sigma_{c}\simeq\dfrac{B^{2}}{4\pi\rho\rm{c^{2}}}.$ (14)

We write the enthalpy of equation (9) of a flow from equations 10 and 13 as

$\omega=\rho{\rm{c^{2}}}(1+\sigma)(1+h).$ (15)

We can plug this equation into equation (5) to calculate the particle number
density at the jet base

$n_{0}=\frac{L_{\rm{\rm jet}}}{2\beta_{0}\gamma_{0}c\pi r_{0}^{2}\,({\rm
m_{p}}/\eta_{e}+{\rm m_{e}})c^{2}(1+\sigma_{\rm c})}.$ (16)

We further use the relativistic Bernoulli’s equation to express the
conservation of energy flux along the jet axis (Königl, 1980)

$\gamma\dfrac{\omega}{\rho}=\rm constant,$ (17)

and from equation (15) we rewrite the above equation such as to define:

$\mu\equiv\gamma\,(1+\sigma)(1+h),$ (18)

where $\mu$ is the normalised total energy flux and is conserved along the
jets (unless the jets entrain mass; see below). In a cold jet where the
specific enthalpy $h$ is negligible, equation (18) simplifies to
$\mu\simeq\gamma\,(1+\sigma_{c})$. This is a very well-known equation to
express the maximum jet Lorentz factor when the majority of the Poynting flux
has been converted to kinetic energy ($\gamma_{\rm max}\simeq\mu$). In this
work, we keep this term in our calculations because $h$ is an estimate of the
energy that the accelerated particles carry in each jet segment, and in
numerous instances can dominate both the magnetisation and the jet Lorentz
factor.

While the jets accelerate between the launching point and the acceleration
region $z_{\rm acc}$, $\mu$ remains constant. We write equation (18) at the
jet base and equate it to the acceleration region and solve for the initial
magnetisation

$\begin{split}&\gamma_{0}(1+\sigma_{0})(1+h_{0})=\gamma_{\rm
acc}(1+\sigma_{\rm acc})(1+h_{\rm acc})\Rightarrow\\\
&\sigma_{0}=\frac{\gamma_{\rm acc}}{\gamma_{0}}\left(1+\sigma_{\rm
acc}\right)\left(\dfrac{1+h_{\rm acc}}{1+h_{0}}\right)-1,\end{split}$ (19)

and in general for every z below the acceleration region

$\sigma(z\leq z_{\rm
acc})=\frac{\gamma_{0}}{\gamma}\left(1+\sigma_{0}\right)\left(\dfrac{1+h_{0}}{1+h}\right)-1,$
(20)

or

$\sigma(z\leq z_{\rm acc})=\frac{\gamma_{\rm acc}}{\gamma}\left(1+\sigma_{\rm
acc}\right)\left(\dfrac{1+h_{\rm acc}}{1+h}\right)-1.$ (21)

With the magnetisation and the specific enthalpy at the acceleration region as
free parameters ($\sigma_{\rm acc}$ and $h_{\rm acc}$, respectively), we set
the initial magnetisation $\sigma_{0}$ required for the flow to be Poynting
flux dominated and to carry enough energy to efficiently accelerate particles
to non-thermal energies. In particular, we use $\sigma_{\rm acc}$ as a free
parameter because this is the simplest way to force our semi-analytical model
to have dissipated the majority of the magnetisation at the acceleration
region, and we set $h_{\rm acc}$ from equation (10) (see also the discussion
on particle acceleration below). The initial specific enthalpy $h_{0}$ is set
by the free parameters at the jet base, and as we discuss below, it is
negligible for the standard case of an initially cold jet that we study here
(see subsection 3.1).

Above the acceleration region, we assume the toroidal component dominates the
poloidal component of the magnetic fields similar to Blandford & Königl
(1979), so

$B(z>z_{\rm acc})=B_{\rm acc}\left(\frac{z}{z_{\rm acc}}\right)^{-1},$ (22)

where $B_{\rm acc}$ is the magnetic field strength at the acceleration region.

Based on equation (13), we generalize the expression of $\sigma$ for every $z$
above the acceleration region

$\sigma(z\geq z_{\rm acc})=\sigma_{\rm acc}\frac{\rho_{\rm acc}(1+h_{\rm
acc})}{\rho(1+h)}\left(\dfrac{z}{z_{\rm acc}}\right)^{-2}.$ (23)

###  2.2 The acceleration region and particle acceleration

We assume that the pairs at the jet base follow a Maxwell-Jüttner distribution
(MJ; the relativistic regime of the Maxwell-Boltzmann distribution) with a
peak energy $k_{\rm B}T_{\rm e}$ that is a free parameter. The population of
protons on the other hand is cold, making the flow cold at the launching
point.

By the time the flow reaches the acceleration region the Poynting flux
dominated flow has dissipated the magnetic energy, hence the magnetisation has
dropped to a value $\sigma_{\rm acc}$. At the same region, we assume a
constant fraction $f_{\rm pl}\sim 0.1$ of particles accelerates to a non-
thermal power law between a minimum and a maximum energy. For the leptonic
scenario, we assume that only pairs accelerate in a power law from an energy
$\varepsilon_{\rm min}{\rm m_{e}c^{2}}=k_{\rm B}T_{\rm e}$ to some
$\varepsilon_{\rm max}$ that we calculate self-consistently by equating the
acceleration timescale $4\varepsilon{\rm m_{e}c^{2}}/(3f_{\rm sc}ecB)$ to the
escape timescale (Jokipii, 1987; Aharonian, 2004). The acceleration efficiency
$f_{\rm sc}$ depends on the particle acceleration mechanism, but we fix it at
a value between 0.01 and 0.1 leading to a maximum electron energy of the order
of GeV for the case of a BHXB. For the lepto-hadronic scenario, we assume that
protons accelerate as well in a power law from an $\varepsilon_{\rm min}$= 1
to some $\varepsilon_{\rm max}$ that we calculate by equating the acceleration
timescale to the (lateral) escape timescale $r$/c of the jet segment and for
the case of BHXBs it may attain values of the order of 100 TeV and above (Pepe
et al., 2015; Kantzas et al., 2021; Kantzas et al., 2022). We constrain the
non-thermal particle distributions by assuming that they extend up to the
maximum energy, and then they drop exponentially

$\dfrac{{\rm{d}}n\left(\varepsilon\right)}{{\rm{d}}\varepsilon}=K\varepsilon^{-p}\,\exp{\left(-\varepsilon/\varepsilon_{\rm{max}}\right)},$
(24)

where $n$ is the particle number density for any species, $K$ is the
normalisation, and the slope $p$ of the power law depends on the particle
acceleration mechanism, but we use it as a free parameter between 1.7 and 2.4,
assuming it remains the same between electrons and protons.

Finally, we derive the average Lorentz factor for every species from the
equation

$\langle\varepsilon\rangle=\dfrac{\bigintss\varepsilon\dfrac{{\rm d}n}{{\rm
d}\varepsilon}{\rm d}\varepsilon}{\bigintss\dfrac{{\rm d}n}{{\rm
d}\varepsilon}{\rm d}\varepsilon}.$ (25)

###  2.3 Jet evolution and particle acceleration

Beyond the acceleration region where particles accelerate to non-thermal
energies as well, the specific enthalpy can become important because the
average Lorentz factors of pairs and/or protons may have significantly
increased (see equation 10). We write the bulk Lorentz factor for every jet
segment above the acceleration region for an outflow from equation (18):

$\gamma(z)=\gamma_{\rm acc}\left(\dfrac{1+h_{\rm
acc}}{1+h}\right)\left(\dfrac{1+\sigma_{\rm acc}}{1+\sigma}\right).$ (26)

Parameter | Units | Fiducial value(s) | Definition | Equation
---|---|---|---|---
$z$ | $r_{g}$ | $-$ | distance from the black hole along the jet axis | $-$
$z_{0}$ | $r_{g}$ | 6 | distance of the jet base from the black hole | $-$
$\gamma$ | $-$ | $1-3$ | bulk Lorentz factor of the flow | 1
$\gamma_{0}$ | $-$ | 1.1 | bulk Lorentz factor at the jet base | $-$
$r$ | $r_{g}$ | $-$ | cross-sectional radius of the flow | 2
$\theta$ | rad | $-$ | jet opening angle | 3
$n$ | $\rm cm^{-3}$ | $-$ | jet (total) particle number density | 4
$n_{0}$ | $\rm cm^{-3}$ | $-$ | jet number density at the jet base | 16
$n_{\rm e}$ | $\rm cm^{-3}$ | $-$ | jet pair number density | $-$
$n_{\rm p}$ | $\rm cm^{-3}$ | $-$ | jet proton number density | $-$
$\rho$ | $\rm g\,cm^{-3}$ | $-$ | jet mass density | 7
$\omega$ | $\rm erg\,cm^{-3}$ | $-$ | total jet enthalpy | 9
$h$ | $-$ | $-$ | jet specific enthalpy | 10
$\sigma$ | $-$ | $-$ | magnetisation of the flow | 13
$\sigma_{0}$ | $-$ | $1-100$ | magnetisation of the flow at the jet base | 19
$\mu$ | $-$ | $1-100$ | normalised total jet energy flux | 18
$\langle\varepsilon_{\rm e,p}\rangle$ | $-$ | $1-100$ | particle average Lorentz factor | 25
$z_{\rm acc}$ | $r_{g}$ | $10^{3}$ | location where jet acceleration reaches the max value | free parameter
$\gamma_{\rm acc}$ | $-$ | 3 | maximum Lorentz factor of the flow at $z_{\rm acc}$ | free parameter
$r_{0}$ | $r_{g}$ | $10-10^{2}$ | jet base radius | free parameter
$L_{\rm jet}$ | $L_{\rm Edd}$ | 0.002-0.02 | injected jet power at the jet base | free parameter
$\eta_{e}$ | $-$ | $1-10^{6}$ | jet pair-to-proton content | free parameter
$\sigma_{\rm acc}$ | $-$ | 0.1 | magnetisation of the flow at the acceleration region | free parameter
$k_{B}T_{e}$ | keV | $-$ | electron peak energy at the jet base | free parameter
Table 1: The definition of the jet quantities we use in this work with their
units, some fiducial values (if applicable), the equation number where we
define the parameter or whether it is a free parameter. See Sections 2 and 4
for further information.

###  2.4 Radiative Processes

We suggest the interested readers to seek for further details on the radiative
processes in Lucchini et al. (2022) for the leptonic processes, and in Kantzas
et al. (2021) for the hadronic processes. We nevertheless briefly discuss the
main processes here for completeness.

####  2.4.1 Leptonic processes

The main three radiative processes of leptonic nature that we require in our
analysis here are: synchrotron radiation, inverse Compton scattering (ICS) and
pair production. In particular, the thermal pairs of the MJ distribution and
the non-thermal power-law tail above the dissipation region, lose energy due
to cyclo-synchrotron radiation (Blumenthal & Gould, 1970; Rybicki & Lightman,
2008). We only account for the average magnetic field strength of the
particular jet segment and assume an isotropic distribution of pitch angles
that we average over.

We further account for the ICS between the pairs and the radiation fields of
the outflow (Blumenthal & Gould, 1970; Rybicki & Lightman, 2008). In
particular, in this work we neglect any external photon field and only allow
for ICS between the emitting pairs and the synchrotron photons (synchrotron
self Compton; SSC). Plausible external photon fields may be important in the
case of AGN jets but for the study-cases as BHXBs we discuss in this work, we
have shown in previous works that the external photon fields are not critical
(see e.g. Lucchini et al. 2021; Kantzas et al. 2021, however, see also
Zdziarski et al. 2014 and Zacharias et al. 2022 for cases where the external
photon fields may be important to explain the $\gamma$-ray spectrum). For
simplicity, we also neglect any accretion disc in the following discussion,
but we do account for it when examining particular sources, following Lucchini
et al. (2022). For the ICS processes, we account for the Klein-Nishina regime
when necessary, and allow for multiple scatterings to better capture the
evolution of the exponential cutoff. This particular process is the most
computationally expensive amongst the leptonic ones, we hence choose to
neglect it when the radiative output becomes $10^{4}$ times smaller than the
synchrotron counterpart for the particular segment.

The final process of leptonic nature we account for is the photon annihilation
to pair production and vice versa (Coppi & Blandford, 1990). These two
processes are usually negligible, so we do not mention them unless we discuss
their impact on the particle population or the spectrum (see, e.g., Connors et
al., 2019).

####  2.4.2 Hadronic Processes

We account for both proton-proton (pp) and proton-photon (p$\gamma$) processes
when accelerated protons interact with the cold protons of the flow and the
jet radiation, respectively. In particular, we use the semi-analytical
parametrisation of Kelner et al. (2006) for the pp interactions, and Kelner &
Aharonian (2008) for the p$\gamma$. The above analysis provides the resulted
distributions of secondary particles (pions that decay into muons, and the
muons decay into neutrinos, pairs and $\gamma$-rays) and hence cannot account
for any synchrotron radiation of muons and/or pions, but for the current
systems we examine, we see that it is not required. We do however consider the
cyclo-synchrotron radiation of secondary pairs due to the presence of the
magnetic field.

In our particular analysis, we find that the synchrotron photons produced by
the primary pairs act as the target for the p$\gamma$ interactions. Based on
this analysis, we can also produce the neutrino counterpart in a self-
consistent manner (Kantzas et al. in prep).

##  3 Results for the steady-state jets

We first present the results of the analysis of the model where we do not
account yet for any mass entrainment. In this flavour of the model, we try to
better understand and constrain the number of leptons in the jets with respect
to the number of protons $\eta_{e}$. We further present the jet dynamical
properties and their corresponding multiwavelength spectra before we compare
them to ones when we account for mass-loading.

###  3.1 Specific enthalpy and particle acceleration

(a) Purely leptonic acceleration with $\varepsilon_{\rm e,min}=1.5$. (b)
Purely leptonic acceleration with $\varepsilon_{\rm e,min}=10$. (c)
Leptohadronic acceleration with $\varepsilon_{\rm e,min}=1.5$. (d)
Leptohadronic acceleration with $\varepsilon_{\rm e,min}=10$.

(e) More efficient hadronic, and leptonic acceleration with $\varepsilon_{\rm
e,min}=1.5$. (f) More efficient hadronic, and leptonic acceleration with
$\varepsilon_{\rm e,min}=10$.

Figure 1: The jet specific enthalpy $h$ as a function of the jet content
$\eta_{e}=n_{\rm e}/n_{\rm p}$. In all plots, we assume a soft non-thermal
power law with $p=2.2$ to derive the average particle Lorentz factors from
equation (25). The color-map corresponds to the average Lorentz factor of
electrons, with lighter colors to indicate more efficient acceleration. In the
left column and for a less efficient electron acceleration, the minimum
Lorentz factor of the pairs is $\varepsilon_{\rm e,min}=1.5$, whereas in the
right column with a more efficient electron acceleration, $\varepsilon_{\rm
e,min}=10$. In the top panels, we assume only leptonic acceleration, in the
middle, we assume non-efficient hadronic acceleration with $\varepsilon_{\rm
p,min}=1$ and $\varepsilon_{\rm p,max}=100$, and in the bottom panels, we
assume efficient hadronic acceleration with $\varepsilon_{\rm p,min}=10$ and
$\varepsilon_{\rm p,max}=10^{7}$. The vertical lines correspond to
$\eta_{e}=\rm m_{p}/m_{e}$. Overall, the specific enthalpy $h$ may attain
values greater than unity and may hence significantly alter the jet
kinematics. See the text for more details.

We can express equation (10) as

$h=\dfrac{\Gamma_{\rm e}(\langle\varepsilon_{\rm e}\rangle-1)+\Gamma_{\rm
p}(\langle\varepsilon_{\rm p}\rangle-1)\dfrac{{\rm
m_{p}/m_{e}}}{\eta_{e}}}{1+\dfrac{{\rm m_{p}/m_{e}}}{\eta_{e}}},$ (27)

where we used equations (7), (12), and $n_{\rm p}=n_{\rm e}/\eta_{e}$.

From the above equation, we see that the specific enthalpy depends merely on
the ratio between pairs and protons. Moreover, we see that $h$ strongly
depends on any mechanism (acceleration or cooling) that would significantly
change the average Lorentz factor of the particles.

In Fig. 1, we plot the specific enthalpy $h$ as a function of the pair-to-
proton ratio $\eta_{e}$ for various values of $\langle\varepsilon_{\rm
e}\rangle$ and $\langle\varepsilon_{\rm p}\rangle$. Both
$\langle\varepsilon_{\rm e}\rangle$ and $\langle\varepsilon_{\rm p}\rangle$
depend on the power law slope of the accelerated particles, as well as the
minimum and the maximum particle energy. We let $\eta_{e}$ to scale between a
few and $10^{6}$ although the latter values are extreme and perhaps not
physically expected. A jet with more protons than leptons ($\eta_{e}<1$) would
be positively charged and hence is unphysical. On the other hand, a very large
number of pairs per protons would be difficult to explain the observed Lorentz
factors on parsec scales (Ghisellini & Tavecchio, 2010).

In the top left plot of Fig. 1 where no protons accelerate at all, and in
particular in the case of approximately equal amount of pairs and protons
($\eta_{e}\sim 1$), we see that the specific enthalpy is significantly smaller
than unity ($h\ll 1$). This is in agreement with the initial setups of GRMHD
simulations where the specific enthalpy is usually neglected (McKinney, 2006;
Komissarov et al., 2007). In the other regime, where the flow is dominated by
pairs ($\eta_{e}\gtrsim 10^{3}$), we see that
$h\sim\Gamma_{e}\langle\varepsilon_{\rm e}\rangle$ (equation 27). In the top
right plot of Fig. 1 where we assume $\varepsilon_{\rm e,min}=10$, we see a
similar evolution of $\eta_{e}$. The main difference is that
$\langle\varepsilon_{\rm e}\rangle$ goes to larger values, hence $h$ goes to
larger values as well. From both plots, we see that for a purely leptonic
flow, the specific enthalpy is not negligible and in fact, it can be as
important as the magnetisation and the kinetic energy in the evolution of the
jets (as discussed below).

In the middle plots of Fig. 1, where protons accelerate in a similar power law
as the accelerated pairs, we see a significantly different evolution of $h$
for different jet content. In particular, in the case where $\varepsilon_{\rm
e,min}=1$ and $\varepsilon_{\rm p,min}=1$ (middle left plot), we see that for
an equal pair-to-proton jet content ($\eta_{e}=1$), $h$ is driven by the
accelerated protons and in fact, $h\sim\Gamma_{p}\langle\varepsilon_{\rm
p}\rangle$ (see equation 27). In the regime of a purely leptonic flow
($\eta_{e}\gg 1$), we see that $h\sim\Gamma_{e}\langle\varepsilon_{\rm
e}\rangle$ and depending on whether $\langle\varepsilon_{\rm
e}\rangle>\langle\varepsilon_{\rm p}\rangle$ or $\langle\varepsilon_{\rm
e}\rangle<\langle\varepsilon_{\rm p}\rangle$, $h$ will increase or decrease,
respectively. In the right-hand-side of the middle panels of Fig. 1, we get
larger values of $\langle\varepsilon_{\rm e}\rangle$ because of the larger
value of $\varepsilon_{\rm e,min}$ (for the particular $p=2.2$), and hence the
specific enthalpy may attain significantly larger values reaching values of
the order of $\Gamma_{e}\langle\varepsilon_{\rm e}\rangle$.

In the bottom plots of Fig. 1 where protons accelerate in a power law from a
$\varepsilon_{\rm p,min}=10$, we see that a flow of $\eta_{e}\sim 1$ has a
significant fraction of energy in the specific enthalpy because
$h\sim\Gamma_{p}\langle\varepsilon_{\rm p}\rangle\sim 90$. In the purely
leptonic regime ($\eta_{e}\gg 1$), we see that $h$ can drop to values smaller
than 10 depending on the average Lorentz factor of the pairs. In the case
where pairs accelerate in a power law from a high energy as $10\,\rm
m_{e}c^{2}$ (right-hand-side plot of lowermost panels of Fig. 1), the energy
content in the specific enthalpy remains significant for both $\eta_{e}\sim 1$
and $\eta_{e}\sim 10^{6}$.

From Fig. 1, we overall see that the specific enthalpy of a flow that
accelerates particles can be important in the evolution of the flow (see also
discussion below). In the case where only pairs accelerate in the jets and for
an equal amount of electrons-to-protons as is commonly assumed in GRMHD (left-
hand-side of the uppermost panels, and in particular in the case of one), we
see that the specific enthalpy is indeed negligible ($h\ll 1$). In any other
case where both pairs and protons accelerate in the jets, and regardless of
the jet content (either pair-dominated or equal pair-to-proton content), the
specific enthalpy of the flow might be of the order of a few-to-tens, and
hence it is important for the evolution of the flow (see also discussion of
CLTM19).

In the Appendix A we discuss the evolution of $h$ in the case of a hard power
law of accelerated particles with $p=1.7$ power law index. Such hard values,
resulting from efficient particle acceleration e.g., in magnetic reconnecting
regions (Sironi et al., 2015; Ball et al., 2018) or relativistic shocks (see
e.g. Böttcher & Baring, 2019), lead to even larger values of $h$ of the order
of thousands. Such large values of $h$ along with large bulk Lorentz factors
as observed in relativistic outflows in AGN and GRBs, would lead to
significantly larger values of total energy flux $\mu$ compared to those in
the literature (Komissarov et al., 2007, 2009; Petropoulou et al., 2022).
Furthermore, equation 18, which has broadly been used to provide an estimate
for the maximum bulk Lorentz factor when the magnetic energy has been
converted into kinetic energy, would not hold anymore and a more careful
treatment where the specific enthalpy is calculated from first principles is
needed.

###  3.2 Total energy flux evolution for steady state jets

In Fig. 2, we plot the evolution of $\mu$ along the jets with the different
components: magnetisation ($\sigma$), bulk Lorentz factor ($\gamma$) and
specific enthalpy ($h$). In the left plots of Fig. 2, we assume a jet content
of equal number of leptons and protons ($\eta_{e}=1$) and in the right plots
we assume a pair-dominated outflow ($\eta_{e}=10000$). In the top panels, we
assume that only leptons accelerate to non-thermal energies, whereas in the
bottom panels, we assume that hadrons accelerate as well in a power law with
the same index.

In the top left panel, where we account only for leptonic acceleration with
$\langle\varepsilon_{\rm e}\rangle=6$, we see that the initial magnetisation
of the outflow converts to bulk kinetic energy whereas the magnetisation drops
to $\sigma_{\rm acc}=0.1$ (a free parameter). The specific enthalpy starts as
negligible at the cold jet base ($h_{0}\ll 10^{-2}$) and remains insignificant
for the jet evolution above the particle acceleration region $z_{\rm acc}$.
This particular regime where the specific enthalpy is insignificant and the
jet composition is one lepton per proton, is the regime considered by most
GRMHD simulations (see also Section 4), and in fact, is the only regime that
BHJet can probe self-consistently so far (Lucchini et al., 2022). With the
current improvement of this work, we can now further explore the jet
kinematics to other regimes where the distribution of the internal energy
density is important in the evolution of the jet dynamics and the
electromagnetic spectrum.

In the top right panel, where we assume a pair-dominated jet ($\eta_{e}\gg 1$)
that accelerates only leptons, we see that the initial magnetisation converts
almost equally to bulk kinetic energy and internal energy ($h$ is now
comparable to $\gamma$). The initial specific enthalpy at the jet base is
larger compared to the previous case and based on equation (18), we see that
also $\mu$ has significantly increased (see also section 3.1).

In the bottom left panel of Fig. 2 where we account for hadronic acceleration
with $\langle\varepsilon_{\rm p}\rangle=4$, we see that the initial
magnetisation dissipates almost equally to kinetic and internal energy. The
initial specific enthalpy is negligible at the cold jet base but when
particles accelerate at the acceleration region, $h$ increases to values
comparable to $\gamma$. Finally, in the bottom right panel where the jet is
pair-dominated, we see that the specific enthalpy at the jet base is of the
order of 1 but still much smaller than the initial magnetisation.

In Fig. 2, according to the approach we follow here, $h$ can overall be
significant for the jet evolution depending on the hadronic acceleration and
the jet content. The former, in particular, strongly depends on the jet
properties, but we cannot capture this non-linear behaviour of the jet
evolution, its effect on the particle acceleration and the consequent feedback
of particle acceleration back to the jet evolution without significantly
increasing the computational cost of the model. However, we can still
investigate the jet properties to gain a better insight on jet physics.

In Appendix B we present a more detailed series of jet evolution for various
jet quantities and different average particle Lorentz factors. Overall, we
find that for many physical scenarios, the specific enthalpy becomes important
for the jet evolution, especially in the case where hadrons accelerate in the
jets as well, and for pair-dominated outflows (see also section 3.1).

(a) Pair-proton outflow with only leptonic acceleration

(b) Pair-dominated outflow with only leptonic acceleration

(c) Pair-proton outflow with lepto-hadronic acceleration

(d) Pair-dominated outflow with lepto-hadronic acceleration

Figure 2: The energy jet components $\gamma$ (the bulk Lorentz factor),
$\sigma$ (the magnetisation), and $h$ (the specific enthalpy) that follow the
relation $\mu=\gamma(1+\sigma)(1+h)$ (equation 18). In all plots we use
$z_{0}=6\,r_{g}$, $z_{\rm acc}=10^{3}\,r_{g}$, $\gamma_{\rm acc}=3$,
$\sigma_{\rm acc}=0.1$ and $\langle\varepsilon_{\rm e}\rangle=6$ (see Table 1
for definitions). We show a pair/proton flow with $\eta_{e}=1$ in the left
column and a pair-dominated flow with $\eta_{e}=10000$ in the right column. In
the top panels, we only account for leptonic acceleration and in the bottom
panels, we consider hadronic acceleration as well with
$\langle\varepsilon_{\rm p}\rangle=4$. The specific enthalpy $h$ leads to
different jet dynamical quantities based on whether hadronic acceleration
takes place and jet content.

(a) Pair-proton outflow with only leptonic acceleration

(b) Pair-dominated outflow with only leptonic acceleration

(c) Pair-proton outflow with lepto-hadronic acceleration

(d) Pair-dominated outflow with lepto-hadronic acceleration

Figure 3: The predicted spectral energy distributions for the four models of
Fig. 2. In the top panels, we only account for leptonic acceleration, and in
the bottom ones, we consider both leptonic and hadronic. In the two left
plots, we assume one proton per electron ($\eta_{e}=1$) and in the right ones
we assume $\eta_{e}=10^{4}$. In all four panels, we use $k_{B}T_{e}=500\,\rm
keV$, and $z_{\rm diss}=1000\,r_{g}$ for a 10 $\rm M_{\odot}$ BHXB at 3 kpc.
We also assume $L_{\rm jet}=2\times 10^{-2}\,L_{\rm Edd}$ for the leptonic
scenarios and $L_{\rm jet}=2\times 10^{-3}\,L_{\rm Edd}$ for the hadronic. The
aforementioned values lead to $\langle\varepsilon_{\rm e}\rangle=5$. We
highlight the contribution of the jet-segments before the dissipation region
(yellow shaded) and that of the jet-segments above the dissipation region
(blue shaded). We show the synchrotron emission with densely-dashed green
line, and the contribution of the ICS with loosely-dashed blue line. Finally,
the green shaded region is the hadronic contribution where we include both
neutral pion decay and the synchrotron radiation of the secondary electrons.

###  3.3 Electromagnetic spectrum of steady state jets

We plot in Fig. 3 the multiwavelength spectra that correspond to the four
different models of Fig. 2. In particular, in the top panels we plot the
purely leptonic scenarios, whereas in the bottom we plot the lepto-hadronic
models. For the left plots, we assume one proton per electron ($\eta_{e}=1$),
whereas on the right plot we examine the extreme case of $\eta_{e}=10^{4}$.

For all four panels, we assume a quite “warm” MJ distribution of leptons with
$k_{B}T_{e}=1000\,\rm keV$, an initial jet-radius of $10\,r_{g}$ in which we
inject some power equal to $10^{-2}\,L_{\rm Edd}$ for the leptonic models, and
$10^{-3}$ for the lepto-hadronic ones. The particle acceleration that happens
at $1000\,r_{g}$ leads to a power-law of particles with an index of 2.2. In
all panels, we show the contribution to the spectrum of the jet segments
before the dissipation region (yellow-shaded) and above (blue-shaded). For the
lepto-hadronic model of the bottom panels, we include the hadronic
contribution as green-shaded. Finally, the densely-dashed line shows the
synchrotron contribution, whereas the loosely-dashed line corresponds to the
ICS.

In the top left panel of Fig. 3, we see the emission from the thermal
electrons dominates in the UV and X-ray bands, whereas the outer jets dominate
in the radio bands via synchrotron radiation, and in the GeV with ICS. In the
case where we assume an increased ratio of pairs (top right panel), for the
same initial conditions we see once more the emission from the thermal pairs
to dominate the UV/X-ray bands but the X-ray luminosity is increased because
the initial pair number density has increased (see equation 16).

In the lepto-hadronic cases of the bottom panels of Fig. 3, we see that the
pair content may significantly affect the SED, and in particular the high-
energy part. For the case of one proton per lepton, we see that the GeV-to-TeV
spectrum first drops exponentially due to the synchrotron emission of the
primary pairs, but later increases due to the hadronic contribution of the
p$\gamma$ interactions. The ICS contribution in this particular case is well
below the hadronic contribution (loosely-dashed line). In the pair-dominated
jet of the right-hand panel, we see that the increased number of pairs leads
to a stronger GeV-to-TeV flux that dominates over the hadronic contribution.

##  4 Mass loaded jets

Parameter | Fiducial value(s) | Definition | Status
---|---|---|---
$\gamma_{\rm 0}$ | 1.11 | bulk Lorentz factor at the jet base $z_{\rm 0}$ | fixed
$\sigma_{0}$ | $10-50$ | magnetisation of the flow at the jet base | free
$k_{B}T_{e}\,/\rm keV$ | 500 | electron peak energy at the jet base | free
$\gamma_{\rm acc}$ | $2-10$ | bulk Lorentz factor at $z_{\rm acc}$ | free
$h_{\rm acc}$ |  $h_{0}$† | jet specific enthalpy at $z_{\rm acc}$ | fixed
$f_{\rho}$ | 10 | jet mass density increase factor | fixed
$z_{\rm diss}/r_{g}$ | $100$ | region where the mass entrainment initiates‡ | free
$z_{\rm load,end}/z_{\rm diss}$ | $100$ | region where the mass entrainment finishes | fixed
Table 2: The fixed and the free (fitted) parameters that drive the mass-
loading jet dynamics. See Section 4 for further information.
†calculated by the temperature of the electrons at the jet base (see equation
10),
‡same as $z_{\rm acc}$.

High-resolution GRMHD simulations of accreting black holes that launch highly
collimated jets suggest that a significant portion of the wind from the
accretion disc might end up in the jet via entrainment. While the jets
accelerate in a dense surrounding medium, they are subject to lateral pressure
from the wind of the accretion disc that results in jet-wind collisions,
causing the jet to wobble. Pinch instabilities form at the jet-wind interface
close to the black hole, almost independently of the initial magnetisation of
the jet, as long as it starts out Poynting flux dominated. These instabilities
dissipate magnetic energy to heat and increases the specific enthalpy of the
jet (see, e.g., Eichler, 1993; Bowman et al., 1996; Spruit et al., 1997;
Begelman, 1998; Giannios & Spruit, 2006; Bromberg & Tchekhovskoy, 2015).

Interestingly, two properties of a collimated jet change at distances $\sim
10^{2}-10^{3}\,r_{g}$: (1) the toroidal component of the magnetic field starts
to dominate over the poloidal component, and (2) the jet speed exceeds the
local fast magnetosonic wave speed, i.e., becomes superfast, the magnetic
analogue of the fluid becoming supersonic CLTM19. Beyond this region, the jet
becomes more susceptible to instabilities forming at the interface between the
flow and the ambient medium. In particular, magnetic pinch instabilities lead
to the formation of eddies that trap matter from the wind and drive it inwards
through the jet-wind interface, allowing for mass entrainment (Mignone et al.,
2013; Gourgouliatos & Komissarov, 2018; Bodo et al., 2021). Without such
eddies, significant mass entrainment into the jet from the external medium may
not be possible due to the jet’s strong magnetic field. Hence, we link the
region where the mass-loading becomes important because of instabilities
explicitly to the region where non-thermal particle acceleration occurs.
Following the results of CLTM19, we connect this region to the first particle
acceleration region of jets as originally proposed by Markoff et al. (2005);
Polko et al. (2014).

In this work, we parametrise the fiducial model B10 of CLTM19 to derive a
semi-analytical formalism that connects the mass loading region to the
particle acceleration region, and study its impact on the emitted
electromagnetic spectrum by studying both the leptonic and the hadronic
processes we discussed above. We consider B10 for our problem because the jet
undergoes strong collimation out to very large scales. Other models explored
in CLTM19 either have too small an accretion disk, such that there is hardly
any lateral pressure from the disk wind. The jets therefore become
uncollimated and thus conical within $1000\,r_{g}$, and hence do not properly
represent the highly collimated, parabolic, large-scale jets we are targeting.
Further, the fast lateral expansion of the uncollimated jet suppresses pinch
instabilities (Moll et al., 2008; Granot et al., 2011; Porth & Komissarov,
2015) and thus exhibits little to no mass-loading (CLTM19).

CLTM19 confirm that the magnetic energy converts to kinetic energy,
accelerating the jets similar to what was found in previous works (McKinney,
2006; Komissarov et al., 2007, 2009). When matter is entrained by the jets,
further magnetic energy is dissipated to heat up the jet, and the inertia of
the entrained gas slows down the jet. The mass entrainment leads to a decrease
of the total (specific) energy flux $\mu$ along the jets up to the distance
where the mass loading stops. Beyond distances of a few $10^{4}\,r_{g}$, the
CLTM19 jet properties have not achieved steady-state as the slow, mass-loaded
jet is still punching through the ambient medium at this point of time in the
GRMHD simulation. Indeed, the simulations suggest that as the jet slowdown due
to massloading suppresses pinch instabilities further along the jet, and
therefore, massloading becomes considerably weaker beyond $10^{4}\,r_{g}$. As
a result, when the simulated jets attain steady-state out to $\gtrsim
10^{5}\,r_{g}$, we expect that $\mu$ should be conserved for the rest of the
jets and there would be jet re-acceleration while both the magnetisation and
the specific enthalpy decrease.

Inspired by the simulation results, our semi-analytical “mass-loaded” jet
model assumes that the mass loading initiates at a distance $z_{\rm diss}\sim
100\,r_{g}$ and ends at $100\,z_{\rm diss}$, with a net increase of
$f_{\rho}=10$ in the jet mass density. Beyond $100\,z_{\rm diss}$, we assume a
constant $\mu$ and steady jet acceleration. In Fig. 4, we plot the mass
density of a mass loaded jet (solid line) and compare it to a non-loaded
steady state jet, assuming one proton per lepton. We show the resulting energy
components ($\gamma$, $\sigma$, and $h$) of the B10 model of CLTM19 in Fig. 5
with dashed lines, and below, we discuss the way we parametrise these
quantities.

Figure 4: The mass density profile of a mass loaded jet (solid line) compared
to a steady state jet without mass loading (dashed line). Both profiles are
normalized to the initial mass density at the jet base. The mass loading
initiates at a distance $z_{\rm diss}$ and at $100\,z_{\rm diss}$ the mass
density has increased by a factor of 10 compared to a non-loading, steady-
state jet.  Figure 5: The energy flux components of a mass loaded jet, where
$\mu$ is the ratio between the total energy flux and the rest-mass flux,
$\gamma$ is the bulk Lorentz factor, $\sigma$ is the magnetisation, and $h$ is
the specific enthalpy. The mass entrainment occurs between $10^{2}$ and
$10^{4}\,r_{g}$ (vertical lines), but the entrained matter becomes comparable
to the mass of the jet at a distance of $10^{3}\,\rm r_{g}$ (middle vertical
line). Finally, we over-plot with dashed lines the fiducial model B10 of
CLTM19 on which we base our analysis (see Section 4).

###  4.1 Mass loading region

In this section, we present the parametrisation of the values of $\sigma$,
$\gamma$ and $h$ of the mass loading region based on the B10 model described
above. In particular, we fit a polynomial to the CLTM19 profiles along the jet
between $z_{\rm diss}=100\,r_{g}$ and $z_{\rm load,\,end}=10^{4}\,r_{g}$. The
profiles of the three quantities $\sigma$, $\gamma$ and $h$ are hard to
predict in such a complex and non-linear system, we hence decide to fit only
for these three quantities and derive $\mu$ from equation (18).

$\begin{split}\log_{10}(\sigma)=0.621\,x^{5}-3.005\,x^{4}+4.599\,x^{3}\\\
-2.502\,x^{2}+0.242\,x+0.563,\end{split}$ (28)
$\begin{split}\log_{10}(\gamma)=-0.276\,x^{6}+1.412\,x^{5}-2.207\,x^{4}+0.853\,x^{3}\\\
+0.257\,x^{2}-0.075\,x+0.394,\end{split}$ (29)
$\begin{split}\log_{10}(h)=0.467\,x^{5}-1.903\,x^{4}+1.100\,x^{3}\\\
+2.482\,x^{2}-1.171\,x-1.826,\end{split}$ (30)

where $x=\log_{10}\left(z/z_{\rm diss}\right)$ and $1\leq
x\leq\log_{10}\left(z_{\rm diss}/z_{\rm load,\,end}\right)$.

We connect the jet base to the mass-loading region assuming that the specific
enthalpy is constant to its initial value at the jet base as we calculate it
with equation (10). We assume that the flow is launched at some speed equal to
the speed of sound (see equation 1) and reaches a value $\gamma_{\rm acc}$,
which is a free parameter, following a logarithmic dependence. In Table 2, we
show the parameters of the mass-loading jet model, indicating whether they are
fixed or fitted parameters.

###  4.2 Jet segments beyond the mass loading region

Given our assumption that once mass-loading stops, the total energy flux is
again conserved, i.e., $\mu$ is constant. Thus we fix $\mu$ at its value at
the end of the mass-loading region, and to better constrain the profile of
$\sigma$ and $h$ beyond the mass-loading region, we fit a first order
polynomial between $10^{4}$ and $10^{5}\,\rm r_{g}$, with coefficients:

$\log_{10}(\sigma)=-0.097\,x-0.178,$ (31) $\log_{10}(h)=-0.245\,x-0.576,$ (32)

where $x$ is the same as above. Here we choose to interpolate the profile of
$\gamma$ and $\sigma$ from the simulation data that closely follows the
expected slow acceleration profile seen in semi-analytical MHD solutions of
particle-dominated ($\sigma\lesssim 1$) jets (see e.g., Tchekhovskoy et al.,
2009).

Having derived the values of $\mu$, $\sigma$ and $h$, we calculate the bulk
Lorentz factor for every jet segment above the $z_{\rm diss}$

$\gamma(z\geq z_{\rm diss})=\frac{\mu}{\sigma+h+1}.$ (33)

###  4.3 Particle acceleration and mass loaded jets

At the location where matter is entrained into the jets, particles start to
accelerate to non-thermal energy as well. Based on the definition of $h$, we
solve for the energy density of the protons

$U_{p}=\frac{h\rho c^{2}-\Gamma_{e}U_{e}}{\Gamma_{p}},$ (34)

where we calculate $U_{e}$ from equation (11) for an MJ+non-thermal power-law
distribution of electrons with a fixed ratio of 10 per cent between the
thermal and the non-thermal electrons, and a fixed power-law slope $p$. We
finally, derive the normalisation of the non-thermal protons

$K_{p}=\frac{U_{p}}{\rm{m_{p}c^{2}}\int\varepsilon^{-p+1}\exp(-\varepsilon/\varepsilon_{\rm
max}){\rm d}\varepsilon},$ (35)

where

$\frac{{\rm d}n_{p}}{{\rm
d}\varepsilon}=K_{p}\varepsilon^{-p}\exp(-\varepsilon/\varepsilon_{\rm max}).$
(36)

Following the above approach, we manage to self-consistently connect the mass-
loading that leads to an increase of the specific enthalpy $h$ to the
electromagnetic radiation due to the proton acceleration.

##  5 results for mass-loaded jets

###  5.1 Total energy flux evolution for mass-loaded jets

In Fig. 6, we present the energy components for two different mass-loaded jets
following the prescription of Section 4. We assume that both jets are Poynting
flux dominated at the jet base with an initial magnetisation of
$\sigma_{0}=40$ and accelerate to a bulk Lorentz factor of $\gamma_{\rm
acc}=3$. In the left panel, we assume one electron per proton at the jet base
$\eta_{e}=1$, and in the right, we assume a pair-dominated jet of
$\eta_{e}=10000$. In both cases, we set the temperature of the thermal
electrons at the jet base at $k_{B}T_{e}=200\,\rm keV$. In the particular case
of the pair-dominated jets, the specific enthalpy reaches values that are
comparable to or even exceeding that of the bulk Lorentz factor and the
magnetisation, especially at the loading region (see also equation 27).
Despite the initially pair-dominated jet base, the matter entrained into the
jets is in approximately equal number of electrons and protons because we
assume that the most likely composition of an accretion disc wind is a neutral
gas of electrons and protons. The jet composition hence changes from pair-
dominated at the regions before the loading to almost equal number of protons
and pairs (Anglés-Castillo et al., 2020).

In the right panel of Fig. 6, we see that the increased number of pairs at the
jet base leads to an increase of $h$. In the extreme case where $\eta_{e}\gg
1000$, the peak of the profile of $h$ may lead to an artificial and unphysical
increase of $\mu$ in the loading region. In Appendix C we discuss how we
constrain the increase of $h$ to avoid such an artificial “mass loss”.

###  5.2 Electromagnetic spectra of steady state mass-loaded jets

In Fig. 7, we plot in the left the predicted SED of the fiducial mass-loaded
jet model based on the dynamical quantities that we show in Fig. 5. We further
assume a jet base of radius $10\,r_{g}$, an electron temperature of $200\,\rm
keV$ at the jet base, an injected jet power of $10^{-3}\,L_{\rm Edd}$ and the
power-law slope of the accelerated particles $p=2.2$ for both leptons and
protons. Similar to above, we show the contribution of the leptonic emission
of the jet segments before the dissipation/loading region, the leptonic
contribution from the dissipation/loading region and beyond, and the hadronic
contribution that is due to p$\gamma$. In the right subplot, we show the
spectrum of a non-loaded jet with similar initial conditions. The main
differences are in the jet emission from the jet base (yellow-shaded region)
and the hadronic contribution. The jet-base emission is higher in the non-
loaded case due to the magnetisation profile we assume here that leads to
greater values for the first few jet segments up to the acceleration region
(see, e.g., Fig. 2).

In Fig. 8 we plot the SEDs that correspond to the two models of Fig. 6, where
we account for mass loading at a distance $100\,r_{g}$ and we assume the
injected jet power to be $10^{-3}\,L_{\rm Edd}$. We show the different
components of the spectrum in Appendix D.

Figure 6: Similar to Fig. 5, but for: left an initial lepton temperature at
the jet base of $k_{B}T_{e}=200\,\rm keV$ and one electron per proton
($\eta_{e}=1$), and in the right for a pair-dominated jet ($\eta_{e}=10000$).
Both scenarios are for an initial magnetisation of $\sigma_{0}=40$ and
$\gamma_{\rm acc}=3$. The increased pair content of the right subplot leads to
an increased initial specific enthalpy of the jets.

Figure 7: Left: The predicted spectral energy distribution of a mass-loaded
jet that corresponds to the dynamical quantities of Fig. 5 for a 10 $\rm
M_{\odot}$ BHXB at 3 kpc. We assume a jet base of 200 keV and radius of 10
$r_{g}$. We show the contribution of the jet-segments before the mass loading
(yellow shaded region), and the contribution of the mass-loaded segments of
both leptonic (blue-shaded) and hadronic (green-shaded). The hadronic
contributes includes both the neutral pion decay and the synchrotron radiation
of the secondary electrons/positrons. Right: Similar to the left, but for a
non-loaded jet with similar initial conditions.

Figure 8: Similar to Fig. 7 but for the mass-loaded jets that correspond to
the dynamical quantities of Fig. 6. The overall spectral distribution can
significantly change under the assumption of a pair-dominated jet base
($\eta_{e}=10000$) in the right plot.

##  6 Discussion

###  6.1 Steady state jets

In the first part of this work, we present the analytical jet model that
includes the specific enthalpy in the jet kinematics and the spatial
evolution.

####  6.1.1 Specific enthalpy, particle acceleration and jet evolution

The specific enthalpy $h$ is a good estimate of whether a jet is cold or hot,
with values of $h\ll 1+\sigma$ to indicate a cold flow, and values of
$h\gtrsim 1+\sigma$ to indicate a hot flow. Astrophysical jets launched by
black holes are overall considered cold and strongly magnetised. The majority
of semi-analytical models that focus on the radiative output rather than the
detailed description of the jets, neglect the specific enthalpy for simplicity
(Markoff et al., 2005; Bosch-Ramon et al., 2006; Vila et al., 2012; Zdziarski
et al., 2014). When particles accelerate though, and in particular in the case
where these accelerated particles carry a significant fraction of the jet
energy, the specific enthalpy increases. As we show in Fig. 1, the exact value
of the specific enthalpy may get values that can easily compare to the bulk
Lorentz factor (values of the order of 1 and $\sim 10$) and/or the jet
magnetisation (values greater than unity for a magnetised outflow). The exact
value of $h$ strongly depends on three aspects: the matter composition of the
jet, the efficiency of the leptonic acceleration, and whether hadrons
accelerate as well or not.

##### Leptonic acceleration

In the case where only leptons accelerate inside the jets, we expect the
electron average Lorentz factor to increase as the acceleration efficiency
increases (top subplots of Fig. 1) and hence the specific enthalpy to increase
as well, according to equation (27). The total specific enthalpy however
depends on the jet composition as well. When a jet is of one electron per
proton ($\eta_{e}=1$), the values of $h$ are $\sim 0.01$ regardless of the
exact average Lorentz factor of the electrons (as long as the average Lorentz
factor of the leptons remains less than $\rm m_{p}/m_{e}\simeq 1836$). This is
the typical scenario that current GRMHD and semi-analytical jet models
consider when studying the exact jet evolution, both in space and time. As we
mention above though, based on observations of both extragalactic and Galactic
jets, it is very likely that jets are pair-dominated (or at least the scenario
of one proton per electron is disfavoured in some cases). Such a jet content
leads to an increase in the specific enthalpy compared to the case of
$\eta_{e}=1$ (see equation 27). The specific enthalpy hence of a jet that is
pair-dominated at launching may contribute significantly to the spatial
evolution of the jet, and the more relativistic (or warmer) the distribution
of pairs, the larger the impact of $h$ on the jet evolution. A pair dominated
jet in fact requires specific enthalpy that can be two to three orders of
magnitude larger than the jet case of an equal number of electrons and protons
(see, e.g., top plots of Fig. 1). Consequently, to achieve bulk flow
acceleration up to the same bulk Lorentz factor, a pair-dominated jet, also
requires a larger value of magnetisation at the jet base if energy flux is
conserved along the jet.

##### Lepto-hadronic acceleration

The energy content of the particles can further increase when jets accelerate
both leptons and hadrons to non-thermal energies. In fact, the more efficient
the particle acceleration, the larger the specific enthalpy, which may get
values of the order of $\Gamma_{p}\langle\varepsilon_{p}\rangle$, regardless
of the jet content, as long as $\eta_{e}\leq\rm m_{p}/m_{e}$ (Fig. 1). It is
hard to predict the exact value of the specific enthalpy in a jet that
efficiently accelerates particles, but overall, it may get values equal to or
even exceed that of the bulk Lorentz factor and/or the magnetisation, that
would mean that the outflow converts to particle dominated instead. We hence
suggest that the specific enthalpy should be treated with extra care and
should not always be considered negligible.

####  6.1.2 Specific enthalpy and spectrum

The SED of the steady jets strongly depends not only on the hadronic
acceleration (or lack of it), but also on the jet content. The most important
difference is in the GeV-to-TeV spectrum. A pair dominated jet is
characterised by the ICS and any contribution from the hadronic processes is
suppressed. In the case of a jet with equal number of protons and pairs, and
accounting for an efficient hadronic acceleration, the hadronic component
dominates in the GeV/TeV bands via the neutral pion decay, which has a
distinguishable shape than that of ICS in the Klein-Nishina regime.

The IR-to-X-ray spectrum of BHXBs may be contaminated by different components,
such as the companion star and/or the accretion disc. In the case of a pair-
dominated jet, though, the X-ray spectrum shows the multiple Compton
scatterings due to the increased electron density that can potentially
replicate the role of the theoretical corona (Markoff et al., 2005, 2015;
Lucchini et al., 2021; Cao et al., 2021). Such an X-ray signature can prove a
useful tool to distinguish between different jet compositions, especially with
the next-generation X-ray telescopes, such as for instance the Imaging X-ray
Polarimetry Explorer (IXPE; Weisskopf et al. 2016), the Advanced Telescope for
High-energy Astrophysics (Athena; Nandra et al. 2013) and the Advanced X-ray
Imaging Satellite (AXIS; Mushotzky et al. 2019).

###  6.2 Mass-loaded jets – HadJet

The initial jet composition at the jet base significantly alters the specific
enthalpy of the jet along its axis, even if we assume that at the mass-loading
region the jet converts to a pair-proton outflow. We see, in particular, that
a pair-proton jet base with a thermal pair distribution that peaks at some
energy of the order of 500 keV, which is a reasonable value for BHXBs,
resulting in insignificant specific enthalpy compared to the rest energy
components, namely the magnetisation and the bulk Lorentz factor (see, e.g.,
the left subplot of Fig. 6). If the jet base, on the other hand, is pair-
dominated, then similar to our discussion above, the initial specific enthalpy
at the jet base is increased and hence its effect on the jet dynamical
evolution might be more important because the energy content carried by
particles might be similar to the bulk kinetic energy (see, e.g, the right
subplot of Fig. 6).

The initial conditions at the jet base have a significant impact on the
electromagnetic spectrum that is our tool to distinguish between the different
scenarios. For the two different scenarios we study here, where the one shows
a pair-proton jet base and the other a pair dominated jet base, there are two
prominent differences in the multiwavelength SEDs. The most important one is
in the GeV/TeV regime, where the larger specific enthalpy of the initially
pair-dominated jet base allows for more energy to be transferred to protons.
The increased energy available for non-thermal proton acceleration allows for
a stronger TeV flux, which is dominated by the neutral pion decay due to
p$\gamma$ interactions. Such TeV flux, depending on the distance of the BHXB
(see, e.g., Kantzas et al., 2022) might be significant to be detected by
current TeV facilities, such as the Large High Altitude Air Shower Observatory
(LHAASO) or future $\gamma$-ray facilities, such as the Cherenkov Telescope
Array (CTA). The fact that an initially pair dominated jet can potentially
lead to a stronger TeV flux may sound counterintuitive, but in fact it is
natural in our treatment due to the assumption that the mass loading is linked
with energy dissipation into particle acceleration. The increase of the
specific enthalpy depends on the initial conditions of the jet launching, and
in this work we base our formalism on one specific setup of GRMHD simulations.
A different setup is very likely to lead to less efficient heating of the
jets, the specific enthalpy nevertheless will still increase due to energy
transfer (see discussion of CLTM19). To explore the full range of possible
physical scenarios with GRMHD simulations is currently too computationally
expensive. We can however examine semi-analytically how the impact on the jet
kinematics depends on the level of dissipation by replacing the heating
parameter $f_{\rm heat}$ (that was used in previous work to estimate the
heating of the thermal particles at the particle acceleration region; see,
e.g., discussion in Lucchini et al. 2021) with the fraction of the magnetic
energy that is additionally allowed to go into energising particles. With such
a parameterisation, $h$ will increase by a factor $f_{\rm heat}\sigma$ along
the jet, whereas the magnetisation will be reduced as $(1-f_{\rm
heat})\sigma$. We show in Fig. 9 the impact of this free parameter in the
energy components. To avoid a steep increase of $h$ that looks like a step-
function, we use a function $\tanh^{2}{(z/z_{\rm diss})}$, instead. The
underlying model is that of the left-hand panel of Fig. 6 where we assume a
“hot” jet base (500 keV) and one proton per electron.

A further spectral difference between a pair-proton and a pair-dominated jet
base is in the lower energy regime of the spectrum, and in particular, in the
UV-to-X-ray spectrum. For the same initial magnetisation and injected power,
the number density of the pairs at the pair-dominated jet base is enhanced
(see, e.g., equation 16) resulting in increased Compton scatterings that lead
to a significant difference in the $\sim$1–$100\,\rm keV$ range. The X-ray
spectrum in particular shows a hard spectral index ($\nu
F_{\nu}\propto\nu^{-\alpha+1}$, with $\alpha<1$; see right-hand plot in Fig.
8) that is similar to the expected output of a thermal corona (Sunyaev &
Titarchuk, 1980; Haardt & Maraschi, 1993; Titarchuk, 1994; Narayan & Yi, 1994;
Magdziarz & Zdziarski, 1995).

Figure 9: Similar to the left sub-plot of Fig. 6, but with different $f_{\rm
heat}$ parameters as shown in the plot. The $f_{\rm heat}$ parameter expresses
the fraction of the magnetic energy that is allocated to the specific enthalpy
to allow a further exploration of dissipation beyond our single GRMHD-based
paramerisation.

###  6.3 Proton energy crisis

With the conserved, mass-loading jet model we develop here, we are able to
constrain the total energy that is allocated to the protons and is used to
accelerate them to non-thermal energies. In that way, the total energy carried
by the accelerated protons never exceeds the available energy of the jets that
has been a major issue in the past (Böttcher et al., 2013; Zdziarski &
Böttcher, 2015; Liodakis & Petropoulou, 2020; Kantzas et al., 2022). In Fig.
10, we plot the specific enthalpy of the protons $\Gamma_{p}U_{p}/\rho c^{2}$
divided by $\mu$ as a function of the total jet enthalpy $h$. This quantity
expresses the fraction of the total energy flux of the jet that is used by the
accelerated protons, and we show its dependence on $h$ for different average
electron Lorentz factors, as indicated by the colormap. Regardless of the
average electron energy $\langle\varepsilon_{e}\rangle$, the protons can
hardly carry more than $\sim$10 per cent of the total energy in the jets
because higher fractions would require specific enthalpy $h$ of the order of a
few or above (upper-right corner of the plot) resulting in strongly magnetised
flows ($\sigma\gtrsim\gamma h$). Moreover, for particular values of
$\langle\varepsilon_{e}\rangle$ (see the blue lines for instance that
correspond to values of the order of 1 to 7), the protons can only be
accelerated at $z_{\rm diss}$ and beyond if the total specific enthalpy $h$ is
greater than some critical value $h>h_{\rm crit}$ where

$h_{\rm crit}=\dfrac{(\langle\varepsilon_{e}\rangle-1)\Gamma_{e}}{1+\dfrac{\rm
m_{p}/\rm m_{e}}{\eta_{e}}},$ (37)

hence the cutoffs for different $\langle\varepsilon_{e}\rangle$ at small
values of $h$. In this particular figure, we use $\eta_{e}=10$, but as we show
in Appendix E for smaller (larger) values of $\eta_{e}$ the only difference is
that the cutoffs are located to smaller (larger) values of $h$.

From Fig. 10, we see that the energy of the accelerated protons never exceeds
that of the jet because the specific enthalpy of the non-thermal protons is
always less than the total normalised energy flux ($\Gamma_{p}U_{p}/\rho
c^{2}<\mu$) and hence never violates the energy budget.

Figure 10: The specific enthalpy of the protons $\Gamma_{p}U_{p}/\rho c^{2}$
divided by $\mu$ shows the total energy that is allocated to protons with
respect to the total available jet energy, as a function of the jet specific
enthalpy $h$. We plot the proton energy density for a number of different
electron energy densities that correspond to different values of
$\langle\varepsilon_{e}\rangle$ as shown in the colormap, and we use
$\eta_{e}=10$.

##  7 Summary and conclusions

Relativistic jets are efficient CR accelerators, but we still do not fully
understand the particle acceleration mechanism. To fully interpret the jet
kinematics, and how they relate to particle acceleration, we need to better
understand how to link the observed spectra emitted by jetted sources over
more than ten orders of magnitude in photon frequency to the jet physical
properties. Currently uncertainties about the composition as well as a lack of
conserved dynamical models have contributed to a degeneracy between leptonic
and lepto-hadronic models.

To break this degeneracy, we have developed a new multi-zone approach that
links the jet composition to the jet dynamics. The total energy flux along the
jet is conserved, while magnetic energy can be dissipated into both kinetic
energy and gas enthalpy via particle acceleration. This new approach makes
clear the key role that the specific enthalpy $h$ can have on the evolution
and exchange of energy along the jet. In particular the enthalpy should be
explicitly taken into account in models where: i) electrons accelerate to
large average energies, ii) protons accelerate in the jets as well, and/or
iii) when the jet is pair-dominated, as suggested for numerous Galactic and
extragalactic jets launched by black holes.

When protons are accelerated into a non-thermal power law, the energy
requirement often exceeds the total energy that can be provided by the jet
and/or the accretion energy onto the black hole, potentially violating energy
conservation. We have developed a new model HadJet based on our earlier lepto-
hadronic work, that now conserves energy and includes a prescription for
proton entrainment. Such a mass loading may in fact inhibit proton
acceleration. By allowing the jets to entrain protons over a range of
distance, as seen to occur in GRMHD simulations via eddies forming at the
jet/accretion disc interface (CLTM19), we demonstrate a new method to avoid
the “hadronic power” problem in a more self-consistent approach. In a future
work, we plan to further explore the impact of mass loading on the
multiwavelength emission of both BHXB jets and AGN jets.

## Acknowledgements

We would like to thank the anonymous reviewer for the thorough commenting that
significantly improved the manuscript. DK and SM are grateful for support by
the Netherlands Organisation for Scientific Research (NWO) VICI grant (no.
639.043.513).

## Data availability

No new data were generated or analysed in support of this research.

## References

  * Aartsen et al. (2018) Aartsen M., et al., 2018, Science, 361, 1378
  * Abbasi et al. (2019) Abbasi R. U., et al., 2019, Phys. Rev. D, 99, 022002
  * Abbasi et al. (2020) Abbasi R. U., et al., 2020, ApJ, 898, L28
  * Aharonian (2004) Aharonian F. A., 2004, Very high energy cosmic gamma radiation: a crucial window on the extreme Universe. World Scientific, doi:10.1142/4657
  * Aloisio et al. (2012) Aloisio R., Berezinsky V., Gazizov A., 2012, Astroparticle Physics, 39-40, 129
  * Anglés-Castillo et al. (2020) Anglés-Castillo A., Perucho M., Martí J. M., Laing R. A., 2020, MNRAS, 500, 1512
  * Ball et al. (2018) Ball D., Sironi L., Özel F., 2018, ApJ, 862, 80
  * Begelman (1998) Begelman M. C., 1998, ApJ, 493, 291
  * Begelman et al. (1984) Begelman M. C., Blandford R. D., Rees M. J., 1984, Rev. Mod. Phys., 56, 255
  * Beskin (2010) Beskin V. S., 2010, Physics-Uspekhi, 53, 1199
  * Beskin & Nokhrina (2006) Beskin V. S., Nokhrina E. E., 2006, MNRAS, 367, 375
  * Blandford & Königl (1979) Blandford R., Königl A., 1979, ApJ, 232, 34
  * Blandford & Payne (1982) Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883
  * Blandford & Rees (1974) Blandford R. D., Rees M. J., 1974, MNRAS, 169, 395
  * Blandford & Znajek (1977) Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433
  * Blumenthal & Gould (1970) Blumenthal G. R., Gould R. J., 1970, Reviews of Modern Physics, 42, 237
  * Bodo et al. (2021) Bodo G., Mamatsashvili G., Rossi P., Mignone A., 2021, MNRAS, 510, 2391
  * Bogovalov & Tsinganos (2005) Bogovalov S., Tsinganos K., 2005, MNRAS, 357, 918
  * Bosch-Ramon et al. (2006) Bosch-Ramon V., Romero G. E., Paredes J. M., 2006, A&A, 447, 263
  * Böttcher & Baring (2019) Böttcher M., Baring M. G., 2019, ApJ, 887, 133
  * Böttcher et al. (2013) Böttcher M., Reimer A., Sweeney K., Prakash A., 2013, ApJ, 768, 54
  * Bowman et al. (1996) Bowman M., Leahy J. p., Komissarov S. S., 1996, MNRAS, 279, 899
  * Broderick & Tchekhovskoy (2015) Broderick A. E., Tchekhovskoy A., 2015, ApJ, 809, 97
  * Bromberg & Tchekhovskoy (2015) Bromberg O., Tchekhovskoy A., 2015, MNRAS, 456, 1739
  * Cao et al. (2021) Cao Z., Lucchini M., Markoff S., Connors R. M. T., Grinberg V., 2021, MNRAS, 509, 2517
  * Celotti & Fabian (1993) Celotti A., Fabian A. C., 1993, MNRAS, 264, 228
  * Chatterjee et al. (2019) Chatterjee K., Liska M., Tchekhovskoy A., Markoff S. B., 2019, MNRAS, 490, 2200
  * Connors et al. (2019) Connors R. M. T., et al., 2019, MNRAS, 485, 3696
  * Cooper et al. (2020) Cooper A. J., Gaggero D., Markoff S., Zhang S., 2020, MNRAS, 493, 3212
  * Coppi & Blandford (1990) Coppi P., Blandford R., 1990, MNRAS, 245, 453
  * Corstanje et al. (2021) Corstanje A., et al., 2021, Phys. Rev. D, 103, 102006
  * Crumley et al. (2017) Crumley P., Ceccobello C., Connors R. M. T., Cavecchi Y., 2017, A&A, 601, A87
  * Díaz Trigo et al. (2013) Díaz Trigo M., Miller-Jones J. C., Migliari S., Broderick J. W., Tzioumis T., 2013, Nature, 504, 260
  * Drenkhahn & Spruit (2002) Drenkhahn G., Spruit H. C., 2002, A&A, 391, 1141
  * Eichler (1993) Eichler D., 1993, ApJ, 419, 111
  * Fabrika (2004) Fabrika S., 2004, Astrophys. Space Phys. Reviews, 12, 1
  * Falcke & Biermann (1995) Falcke H., Biermann P. L., 1995, A&A, 293, 665
  * Falcke et al. (2004) Falcke H., Körding E., Markoff S., 2004, A&A, 414, 895
  * Fender et al. (2005) Fender R. P., Maccarone T. J., van Kesteren Z., 2005, MNRAS, 360, 1085
  * Gallo et al. (2005) Gallo E., Fender R., Kaiser C., Russell D., Morganti R., Oosterloo T., Heinz S., 2005, Nature, 436, 819
  * Ghisellini & Tavecchio (2010) Ghisellini G., Tavecchio F., 2010, MNRAS: Letters, 409, L79
  * Ghisellini et al. (2014) Ghisellini G., Tavecchio F., Maraschi L., Celotti A., Sbarrato T., 2014, Nature, 515, 376
  * Giannios & Spruit (2006) Giannios D., Spruit H. C., 2006, A&A, 450, 887
  * Giannios et al. (2004) Giannios D., Kylafis N. D., Psaltis D., 2004, A&A, 425, 163
  * Gourgouliatos & Komissarov (2018) Gourgouliatos K. N., Komissarov S. S., 2018, Nature Astronomy, 2, 167
  * Granot et al. (2011) Granot J., Komissarov S. S., Spitkovsky A., 2011, MNRAS, 411, 1323
  * Guo et al. (2014) Guo F., Li H., Daughton W., Liu Y.-H., 2014, Phys. Rev. Lett., 113, 155005
  * Haardt & Maraschi (1993) Haardt F., Maraschi L., 1993, ApJ, 413, 507
  * Hailey et al. (2018) Hailey C. J., Mori K., Bauer F. E., Berkowitz M. E., Hong J., Hord B. J., 2018, Nature, 556, 70
  * Heinz (2006) Heinz S., 2006, ApJ, 636, 316
  * Heinz & Sunyaev (2003) Heinz S., Sunyaev R. A., 2003, MNRAS, 343, L59
  * Hillas (1984) Hillas A. M., 1984, Annual review of A&A, 22, 425
  * Hjellming & Johnston (1988) Hjellming R., Johnston K., 1988, ApJ, 328, 600
  * Jokipii (1987) Jokipii J., 1987, ApJ, 313, 842
  * Kantzas et al. (2021) Kantzas D., et al., 2021, MNRAS, 500, 2112
  * Kantzas et al. (2022) Kantzas D., Markoff S., Lucchini M., Ceccobello C., Grinberg V., Connors R. M. T., Uttley P., 2022, MNRAS, 510, 5187
  * Kelner & Aharonian (2008) Kelner S., Aharonian F., 2008, Physical Review D, 78, 034013
  * Kelner et al. (2006) Kelner S., Aharonian F. A., Bugayov V., 2006, Physical Review D, 74, 034018
  * Komissarov et al. (2007) Komissarov S. S., Barkov M. V., Vlahakis N., Königl A., 2007, MNRAS, 380, 51
  * Komissarov et al. (2009) Komissarov S. S., Vlahakis N., Königl A., Barkov M. V., 2009, MNRAS, 394, 1182
  * Königl (1980) Königl A., 1980, The Physics of Fluids, 23, 1083
  * Liodakis & Petropoulou (2020) Liodakis I., Petropoulou M., 2020, ApJ, 893, L20
  * Liodakis et al. (2021) Liodakis I., Blinov D., Potter S. B., Rieger F. M., 2021, MNRAS: Letters, 509, L21
  * Lister et al. (2009) Lister M. L., Homan D. C., Kadler M., Kellermann K. I., Kovalev Y. Y., Ros E., Savolainen T., Zensus J. A., 2009, ApJ, 696, L22
  * Lister et al. (2016) Lister M. L., et al., 2016, AJ, 152, 12
  * Lucchini et al. (2018) Lucchini M., Markoff S., Crumley P., Krauß F., Connors R. M. T., 2018, MNRAS, 482, 4798
  * Lucchini et al. (2021) Lucchini M., Russell T. D., Markoff S. B., Vincentelli F., Gardenier D., Ceccobello C., Uttley P., 2021, MNRAS, 501, 5910
  * Lucchini et al. (2022) Lucchini M., et al., 2022, Monthly Notices of the Royal Astronomical Society, 517, 5853
  * Magdziarz & Zdziarski (1995) Magdziarz P., Zdziarski A. A., 1995, MNRAS, 273, 837
  * Maitra et al. (2009) Maitra D., Markoff S., Brocksopp C., Noble M., Nowak M., Wilms J., 2009, MNRAS, 398, 1638
  * Mannheim (1993) Mannheim K., 1993, A&A, 269, 67
  * Mannheim & Schlickeiser (1994) Mannheim K., Schlickeiser R., 1994, A&A, 286, 983
  * Markoff et al. (2001) Markoff S., Falcke H., Fender R., 2001, A&A, 372, L25
  * Markoff et al. (2003) Markoff S., Nowak M., Corbel S., Fender R., Falcke H., 2003, A&A, 397, 645
  * Markoff et al. (2005) Markoff S., Nowak M. A., Wilms J., 2005, ApJ, 635, 1203
  * Markoff et al. (2015) Markoff S., et al., 2015, ApJ, 812, L25
  * Matthews et al. (2020) Matthews J. H., Bell A. R., Blundell K. M., 2020, New Astronomy Reviews, 89, 101543
  * McKinney (2006) McKinney J. C., 2006, MNRAS, 368, 1561
  * Merloni et al. (2003) Merloni A., Heinz S., Di Matteo T., 2003, MNRAS, 345, 1057
  * Mignone et al. (2013) Mignone A., Striani E., Tavani M., Ferrari A., 2013, MNRAS, 436, 1102
  * Mirabel & Rodriguez (1994) Mirabel I., Rodriguez L., 1994, Nature, 371, 46
  * Moll et al. (2008) Moll R., Spruit H. C., Obergaulinger M., 2008, A&A, 492, 621
  * Mori et al. (2021) Mori K., et al., 2021, ApJ, 921, 148
  * Mushotzky et al. (2019) Mushotzky R., et al., 2019, in Bulletin of the American Astronomical Society. p. 107 (arXiv:1903.04083)
  * Nandra et al. (2013) Nandra K., et al., 2013, arXiv:1306.2307
  * Narayan & Yi (1994) Narayan R., Yi I., 1994, arXiv:9403052,
  * Nokhrina et al. (2015) Nokhrina E. E., Beskin V. S., Kovalev Y. Y., Zheltoukhov A. A., 2015, MNRAS, 447, 2726
  * Olejak et al. (2020) Olejak A., Belczynski K., Bulik T., Sobolewska M., 2020, A&A, 638, A94
  * Parfrey et al. (2019) Parfrey K., Philippov A., Cerutti B., 2019, Phys. Rev. Lett., 122, 035101
  * Pepe et al. (2015) Pepe C., Vila G. S., Romero G. E., 2015, A&A, 584, A95
  * Petropoulou et al. (2019) Petropoulou M., Sironi L., Spitkovsky A., Giannios D., 2019, ApJ, 880, 37
  * Petropoulou et al. (2022) Petropoulou M., Psarras F., Giannios D., 2022, Monthly Notices of the Royal Astronomical Society
  * Polko et al. (2014) Polko P., Meier D. L., Markoff S., 2014, MNRAS, 438, 959
  * Porth & Komissarov (2015) Porth O., Komissarov S. S., 2015, MNRAS, 452, 1089
  * Potter & Cotter (2012) Potter W. J., Cotter G., 2012, MNRAS, 429, 1189
  * Protheroe & Kazanas (1983) Protheroe R. J., Kazanas D., 1983, ApJ, 265, 620
  * Pushkarev et al. (2009) Pushkarev A. B., Kovalev Y. Y., Lister M. L., Savolainen T., 2009, A&A, 507, L33
  * Pushkarev et al. (2017) Pushkarev A. B., Kovalev Y. Y., Lister M. L., Savolainen T., 2017, MNRAS, 468, 4992
  * Rachen & Biermann (1993) Rachen J. P., Biermann P. L., 1993, A&A, 272, 161
  * Rachen & Mészáros (1998) Rachen J. P., Mészáros P., 1998, Phys. Rev. D, 58, 123005
  * Reig et al. (2003) Reig P., Kylafis N. D., Giannios D., 2003, A&A, 403, L15
  * Romero et al. (2003) Romero G. E., Torres D. F., Bernadó M. K., Mirabel I., 2003, A&A, 410, L1
  * Rybicki & Lightman (2008) Rybicki G. B., Lightman A. P., 2008, Radiative processes in astrophysics. John Wiley & Sons, doi:10.1002/9783527618170
  * Shikaze et al. (2007) Shikaze Y., et al., 2007, Astroparticle Physics, 28, 154
  * Sironi & Spitkovsky (2014) Sironi L., Spitkovsky A., 2014, ApJ, 783, L21
  * Sironi et al. (2015) Sironi L., Petropoulou M., Giannios D., 2015, MNRAS, 450, 183
  * Sironi et al. (2021) Sironi L., Rowan M. E., Narayan R., 2021, ApJ, 907, L44
  * Spruit et al. (1997) Spruit H. C., Foglizzo T., Stehle R., 1997, MNRAS, 288, 333
  * Sunyaev & Titarchuk (1980) Sunyaev R. A., Titarchuk L. G., 1980, A&A, 86, 121
  * Tchekhovskoy et al. (2008) Tchekhovskoy A., McKinney J. C., Narayan R., 2008, MNRAS, 388, 551
  * Tchekhovskoy et al. (2009) Tchekhovskoy A., McKinney J. C., Narayan R., 2009, ApJ, 699, 1789
  * Tetarenko et al. (2016) Tetarenko B. E., Sivakoff G. R., Heinke C. O., Gladstone J. C., 2016, ApJ Supplement Series, 222, 15
  * Tetarenko et al. (2021) Tetarenko A. J., et al., 2021, MNRAS, 504, 3862
  * The Pierre Auger Observatory et al. (2017) The Pierre Auger Observatory et al., 2017, Science, 357, 1266
  * Titarchuk (1994) Titarchuk L., 1994, ApJ, 434, 570
  * Vila & Romero (2010) Vila G. S., Romero G. E., 2010, MNRAS, 403, 1457
  * Vila et al. (2012) Vila Romero, G. E. Casco, N. A. 2012, A&A, 538, A97
  * Vlahakis & Konigl (2003) Vlahakis N., Konigl A., 2003, ApJ, 596, 1080
  * Waggett et al. (1977) Waggett P., Warner P., Baldwin J., 1977, MNRAS, 181, 465
  * Weisskopf et al. (2016) Weisskopf M. C., et al., 2016, in den Herder J.-W. A., Takahashi T., Bautz M., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 9905, Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray. p. 990517, doi:10.1117/12.2235240
  * Yushkov et al. (2019) Yushkov A., Collaboration P. A., et al., 2019, ICRC, 358, 482
  * Zacharias et al. (2022) Zacharias M., Reimer A., Boisson C., Zech A., 2022, Monthly Notices of the Royal Astronomical Society, 512, 3948
  * Zanin et al. (2016) Zanin R., Fernández-Barral A., de Oña Wilhelmi E., Aharonian F., Blanch O., Bosch-Ramon V., Galindo D., 2016, A&A, 596, A55
  * Zdziarski & Böttcher (2015) Zdziarski A. A., Böttcher M., 2015, MNRAS: Letters, 450, L21
  * Zdziarski et al. (2014) Zdziarski A. A., Pjanka P., Sikora M., Stawarz Ł., 2014, MNRAS, 442, 3243
  * Zdziarski et al. (2022a) Zdziarski A. A., Tetarenko A. J., Sikora M., 2022a, ApJ, 925, 189
  * Zdziarski et al. (2022b) Zdziarski A. A., Phuravhathu D. G., Sikora M., Böttcher M., Chibueze J. O., 2022b, ApJ Letters, 928, L9

##  Appendix A Specific enthalpy for a hard power law of accelerated particles

In Fig. 11 we plot the evolution of $h$ for different jet composition. See
Section 3.1 for a detailed description of the subplots. In this figure, we
show the evolution of $h$ assuming that the particles accelerate in a harder
power law with an index of $p=1.7$ compared to $p=2.2$ we discuss in the main
text.

In the top subplots, we notice a similar behaviour to Fig. 1 but $h$ goes to
larger values for the case of a pair-dominated jet
($h\sim\Gamma_{e}\langle\varepsilon_{\rm e}\rangle$) according to equation
27).

In the case where protons accelerate as well, $h$ can attain values as large
as $\sim 2000$ for a particle acceleration with $\varepsilon_{\rm
e,min}=\varepsilon_{\rm p,min}=10$ as we show in the lowermost subplots. This
value is significantly larger than the expected values of $\gamma$ of the bulk
flow and in combination with the case where $\sigma$ takes large values to
lead to hard power laws of particles (Sironi et al., 2015, 2021; Ball et al.,
2018), we see that the equation $\mu=\gamma(\sigma+h+1)$ (equation 18) would
not be a good approximation for the bulk Lorentz factor anymore (McKinney,
2006; Komissarov et al., 2007, 2009; Beskin, 2010).

(a) Purely leptonic acceleration with $\varepsilon_{\rm e,min}=1.5$. (b)
Purely leptonic acceleration with $\varepsilon_{\rm e,min}=10$. (c)
Leptohadronic acceleration with $\varepsilon_{\rm e,min}=1.5$. (d)
Leptohadronic acceleration with $\varepsilon_{\rm e,min}=10$.

(e) More efficient hadronic, and leptonic acceleration with $\varepsilon_{\rm
e,min}=1.5$. (f) More efficient hadronic, and leptonic acceleration with
$\varepsilon_{\rm e,min}=10$.

Figure 11: The jet specific enthalpy $h$ as a function of the jet content
$\eta_{e}=n_{\rm e}/n_{\rm p}$. In all plots, we assume $p=1.7$ to derive the
average particle Lorentz factors from equation (25). The color-map corresponds
to the average Lorentz factor of pairs with lighter colors to indicate larger
values. In the left column, we set the minimum Lorentz factor of pairs to be
$\varepsilon_{\rm e,min}=1.5$, and on the right column, we use
$\varepsilon_{\rm e,min}=10$. In the top raw, we assume that only leptonic
acceleration takes place, in the middle raw we assume hadronic acceleration as
well with $\varepsilon_{\rm p,min}=1$ and $\varepsilon_{\rm p,max}=100$, and
in the bottom raw, we assume $\varepsilon_{\rm p,min}=10$ and
$\varepsilon_{\rm p,max}=10^{7}$.

##  Appendix B All energy components plots

In Figures 12-15, we show the evolution of $\mu$ along the jet for different
values of h. In particular, for all subplots of Figures 12-15, we assume that
the accelerated particles follow a power law with an index of $p$=2.2. The
outflow launches at a distance of 6 $r_{g}$ from the black hole and the
particle acceleration initiates at $10^{3}\,r_{g}$. While the jets accelerate
at some maximum Lorentz factor $\varepsilon_{\rm acc}=3$, we assume that the
magnetisation at this region has dropped to $\sigma_{\rm acc}=0.1$. For every
subplot, we assume $\eta_{e}=10$ (top left), $\eta_{e}=100$ (top right),
$\eta_{e}=10^{3}$ (bottom left) and $\eta_{e}=10^{4}$ (bottom right) constant
along the outflow.

In Fig. 12, we plot the jet evolution assuming only leptonic acceleration with
an average Lorentz factor of $\langle\varepsilon_{\rm e}\rangle=6$. in
agreement with Fig. 1, we see that while the pair content increases in the
jets, the specific enthalpy increases accordingly, and hence the total $\mu$
increases. In the cases of $\eta_{e}=10^{3}$ and $10^{4}$, in particular, we
see that the specific enthalpy h has values comparable or even larger than the
bulk Lorentz factor of the jet flow.

In Fig. 13, we plot a purely leptonic acceleration similar to Fig. 12 but
assuming $\langle\varepsilon_{\rm e}\rangle=32$ instead. The pair-dominated
jets where $\eta_{3}=10^{3}$ and $10^{4}$ (bottom subplots), indicate that an
efficient acceleration mechanism would lead to high values of h, which for the
case of $\eta_{e}=10^{4}$ the overall value of $\mu$ is of the order of 100, a
much higher value than commonly found in the literature.

In Figures 14 and 15, we further account for hadronic acceleration with
$\langle\varepsilon_{\rm p}\rangle=4$. In the cases where the jets are pair-
dominated, to obtain the specific enthalpy $h$ calculated at the particle
acceleration region, we require a jet base that is Poynting flux dominated
with a magnetisation of the order of 50-100.

(a) $\eta_{e}=10$ (b) $\eta_{e}=100$ (c) $\eta_{e}=1000$ (d) $\eta_{e}=10000$

Figure 12: The evolution of the different energy components $\gamma$, $\sigma$
and $h$ as indicated in each subplot, and the total $\mu$ based on equation
(18). All subplots are for an average electron Lorentz factor of
$\langle\varepsilon_{\rm e}\rangle=6$ and the jet content is shown in each
subplot.

(a) $\eta_{e}=10$ (b) $\eta_{e}=100$ (c) $\eta_{e}=1000$ (d) $\eta_{e}=10000$

Figure 13: Same as Fig. 12 but for $\langle\varepsilon_{\rm e}\rangle=32$.

(a) $\eta_{e}=10$ (b) $\eta_{e}=100$ (c) $\eta_{e}=1000$ (d) $\eta_{e}=10000$

Figure 14: Same as Fig. 12 but for $\langle\varepsilon_{\rm e}\rangle=6$. We
further account for hadronic acceleration with $\langle\varepsilon_{\rm
p}\rangle=4$.

(a) $\eta_{e}=10$ (b) $\eta_{e}=100$ (c) $\eta_{e}=1000$ (d) $\eta_{e}=10000$

Figure 15: Same as Fig. 14 but for $\langle\varepsilon_{\rm e}\rangle=32$ and
$\langle\varepsilon_{\rm p}\rangle=4$.

##  Appendix C Artificial mass loss

In Fig. 16 we show how the increase of the specific enthalpy $h$ in the mass-
loading region may lead to a unphysical increase of $\mu$ that would mean mass
loss instead. Such an artificial mass loss is due to the fact that we assume a
hot flow and/or a pair dominated jet base with $\eta_{e}\gg 1000$. To avoid
such a condition, we first calculate the value of $\mu$ from the 5th order
polynomial

$\begin{split}\log_{10}(\mu)=0.2231\,x^{5}-0.7242\,x^{4}+0.4546\,x^{3}+0.104\,x^{2}\\\
-0.09267\,x+1.031,\end{split}$ (38)

based on the results of CLTM19, and then we calculate $h$ from the equation
$h=\mu/\gamma-(\sigma+1)$, where the values of $\gamma$ and $\sigma$ are from
equations 29 and 28, respectively. In the cases that the specific enthalpy $h$
would drop to zero (see for instance Fig. 16 and 17, we choose to set $h$ to a
very low value, namely one 100th of the magnetisation. In the above equation,
$x$ is the same as in Section 4.

Figure 16: The jet energy components similar to Fig. 6 but for the case of a
jet base with $\eta_{e}=10^{5}$. Following the description we discuss in
Section 4, the particular profile of $h$ leads to an artificial increase of
$\mu$ that would mean mass loss instead, which is unphysical (left). Using the
profile of $\mu$ from equation 38, we constrain $h$ to follow the mass-loading
scenario (right). The initial magnetisation is $\sigma_{0}=5$ and the Lorentz
factor at the dissipation region is $\gamma_{\rm acc}=3$.

Finally, in Fig. 17 we plot the scenario where the artificial mass loss is due
to a combination of a large Lorentz factor ($\gamma_{\rm acc}=10$) and the
profile of $h$ from equation 30. The above assumption allows forcing a mass-
loading scenario.

Figure 17: Similar to Fig. 16 but for the case of $\gamma_{\rm acc}=10$,
$\sigma_{0}=10$ and $\eta_{e}=10^{4}$.

##  Appendix D SED components

In Fig. 18 we show the spectrum as presented on the right subplot of Fig. 8
but with the individual components instead.

Figure 18: Identical to the right subplot of Fig. 8 but in this plot we show
the individual radiative components as indicated in the legend.

##  Appendix E Proton power

In Fig. 19, we plot the fraction of the energy that is allocated to proton
acceleration with respect to the total available energy flux of the jet $\mu$.
We plot this quantity versus the total specific enthalpy of the jet $h$ for
different average Lorentz factors $\langle\varepsilon_{e}\rangle$ of the
electrons. In the main text, we included the case where $\eta_{e}=10$ and here
we plot the cases where $\eta_{e}=1$ (left) and $\eta_{e}=100$ (right), for
completeness. See Section 6.3 for further information.

Figure 19:  The specific enthalpy of the protons $\Gamma_{p}U_{p}/\rho c^{2}$
divided by $\mu$ shows the total energy that is allocated to protons with
respect to the total available jet energy, as a function of the jet specific
enthalpy $h$. We plot the proton energy density for a number of different
electron energy densities that correspond to different values of
$\langle\varepsilon_{e}\rangle$ as shown in the colormap, and we use
$\eta_{e}=1$ in the left, and $\eta_{e}=100$ in the right.