# Modeling particle transport in astrophysical outflows and simulations of
associated emissions

D. A. Papadopoulos,1,3 O. T. Kosmas2 and S. Ganatsios3
1Department of Electrical and Computer Engineering, University of Western
Macedonia, Kozani, Greece
2Modelling and Simulation Center, MACE, University of Manchester, Sackville
Street, Manchester, UK
3Theoretical Physics Section, University of Ioannina, GR-45110 Ioannina,
Greece E-mail<EMAIL_ADDRESS>(DAP)E-mail:
<EMAIL_ADDRESS>(OTK)E-mail<EMAIL_ADDRESS>(SG)

(Last updated 2020 June 10; in original form 2013 September 5)

###### Abstract

In this work, after making an attempt to improve the formulation of the model
on particle transport within astrophysical plasma outflows and constructing
the appropriate algorithms, we test the reliability and effectiveness of our
method through numerical simulations on well-studied Galactic microquasars as
the SS 433 and the Cyg X-1 systems. Then, we concentrate on predictions of the
associated emissions, focusing on detectable high energy neutrinos and
$\gamma$-rays originated from the extra-galactic M33 X-7 system, which is a
recently discovered X-ray binary located in the neighboring galaxy Messier 33
and has not yet been modeled in detail. The particle and radiation energy
distributions, produced from magnetized astrophysical jets in the context of
our method, are assumed to originate from decay and scattering processes
taking place among the secondary particles created when hot (relativistic)
protons of the jet scatter on thermal (cold) ones (p-p interaction mechanism
inside the jet). These distributions are computed by solving the system of
coupled integro-differential transport equations of multi-particle processes
(reactions chain) following the inelastic proton-proton (p-p) collisions. For
the detection of such high energy neutrinos as well as multi-wavelength
(radio, X-ray and gamma-ray) emissions, extremely sensitive detection
instruments are in operation or have been designed like the CTA, IceCube,
ANTARES, KM3NeT, IceCube-Gen-2, and other space telescopes.

###### keywords:

radiation mechanisms: general – neutrinos – binaries: general – stars: winds,
outflows – ISM: jets and outflows – gamma-rays: general

††pubyear: 2020††pagerange: Modeling particle transport in astrophysical
outflows and simulations of associated emissions–A.3

##  1 Introduction

During the last few decades, collimated astrophysical outflows have been
observed to emerge from a wide variety of Galactic and extragalactic compact
structures (Spencer, 1984; Mirabel & Rodríguez, 1999; Romero et al., 2003;
Reynoso, Romero & Christiansen, 2008; Reynoso & Romero, 2009; Romero et al.,
2016). Among such objects, the stellar scale class of microquasar (MQ) and
X-Ray binary systems (XRBs) posses prominent positions (Romero et al., 2016;
Smponias & Kosmas, 2011, 2014; Romero & Vila, 2008; Reid et al., 2011). These
two-body cosmic structures consist of a collapsed stellar remnant, a stellar
mass black hole or a neutron star (compact object), and a companion (donor)
main sequence star in coupled orbit around their center of mass.

Due to the strong gravitational field pertaining around the compact object,
mass from the companion star is accreted onto the equatorial region of the
black hole forming an accretion disc. In the cases when the black hole is
rotating rapidly and the accretion disc is geometrically rather thick and hot,
ejection of two powerful oppositely directed mostly relativistic mass outflows
(jets) occur perpendicular to the accretion disk (Vieyro & Romero, 2017). Such
systems constitute excellent "laboratories" for the investigation of
astrophysical outflows (jets). Currently, for the detection of multi-
wavelength (radio, X-ray and gamma-ray) as well as high energy neutrino
emissions, appreciably sensitive detection tools are in operation or have been
designed like the IceCube, ANTARES, KM3NeT, CTA, IceCube-Gen2, etc. (Saito et
al., 2009; Actis et al., 2011; Aartsen et al., 2015, 2016, 2018; Adrian-
Martinez et al., 2016; Albert et al., 2020).

From a theoretical and phenomenological view point, due to the presence of
rather strong magnetic fields, the MQ jets are treated as magneto-
hydrodymamical flows emanating from the vicinity of the compact object
(usually a stellar mass black hole) (Smponias & Kosmas, 2015, 2017; Kosmas &
Smponias, 2018), and this is also assumed to be the case in our present work.
From the observed characteristics of MQs, researchers have concluded that they
share a lot of similarities in their physical properties with the class of
Active Galactic Nuclei (AGN), which are mostly located at the central region
of galaxies. The latter, however, are enormously larger in scale compared to
microquasars and, in addition, the evolution of AGNs is appreciably slower
than that of MQs that makes the observation of many phenomena significantly
difficult (Remillard & McClintock, 2006; Charles & Coe, 2006; Orosz, 2003;
Fabrika, 2004; Cherepashchuk et al., 2005).

In this article, we focus on the magnetized astrophysical outflows (jets) that
are characterized by hadronic content (p, $\pi^{\pm}$, light nuclei, etc.) in
their jets. We consider them as fluid flows emanating from a central source at
the jet’s origin (Mirabel & Rodríguez, 1999; Smponias & Kosmas, 2011, 2014).
We, furthermore, assume that the compact object is a spinning black hole with
mass up to few tens (30-50) of the Sun’s mass (Smponias & Kosmas, 2014; Romero
et al., 2016; Reynoso, Romero & Christiansen, 2008).

Initially, we make an effort to improve the model employed to perform
numerical simulations for particle and radiation emissions from hadronic MQ
jets (Papadopoulos, 2020; Papadopoulos, Papavasileiou & Kosmas, 2020). We
adopt that the main mechanism producing high energy neutrinos and
$\gamma$-rays, is the proton-proton (p-p) collision taking place within the
hadronic jets, i.e. the inelastic scattering of hot (non-thermal) protons on
thermal (cold) ones (Reynoso, Romero & Christiansen, 2008).

The scattering, diffusion, decay, etc., of the secondary particles
($\pi^{\pm}$, $K^{\pm}$, $\mu^{\pm}$, $e^{\pm}$, etc.) produced afterwards,
from a mathematical modeling point of view, are governed by a system of
coupled integro-differential transport equations. One of the purposes of this
work is to attempt, for a first time, to formulate in a compact way this
differential system of transport equations (satisfied by the primary protons
and the secondary multi-particles, multi-species) and derive advantageous
algorithms to perform the required simulations for the aforementioned
emissivities (Tsoulos, Kosmas & Stavrou, 2017).

Moreover, we are intended to extend the calculations of Ref. (Smponias &
Kosmas, 2014, 2015, 2017; Kosmas & Smponias, 2018) so as to include
contributions to neutrino emissivity (intensity) originating from the
secondary muons ($\mu^{\pm}$) produced from the decay of charged pions
($\pi^{\pm}$) which in our previous works had been ignored due to the long
time consuming required for such simulations. Obviously, the mathematical
problem of such a study becomes also more complicated (e.g. the number of
coupled equations in the above mentioned differential system increases rapidly
in the case of considering neutrinos coming out of kaon ($K^{\pm}$) decays
(Lipari, Lusignoli & Meloni, 2009).

After fixing the model parameters and testing the derived algorithms on the
reproducibility of some known properties of the well studied SS 433 Galactic
micro-quasar (Smponias & Kosmas, 2015, 2017; Kosmas & Smponias, 2018), we
perform detailed simulations for the Galactic Cyg X-1 system as well as for
the extragalactic M33 X-7 MQ (Pietsch et al., 2006). The latter system has not
been studied up to now from a jet emission view point because this is a rather
recently discovered X-ray binary system located in the neighboring galaxy of
our Milky Way Galaxy, known as Messier 33 (Pietsch et al., 2006).

The rest of the paper is organized as follows. In Sect. 2, we describe briefly
the main characteristics of hadronic jets in Galactic micro-quasars (SS433 and
Cyg X-1) as well as the extragalactic M33 X-7 system. Then (Sect. 3), we
present the formulation of our improved method related to the differential
system of transport equations (integro-differential system of coupled
equations) and derive the appropriate algorithms. In Sect. 4 we present and
discuss our results referred to high energy neutrino and $\gamma$-ray energy-
spectra emitted from the extragalactic M33 X-7 MQ. Finally (Sect. 5), we
summarize the main conclusions extracted from the present investigation.

##  2 Brief description of the hadronic MQ systems

In this section we summarize briefly some basic properties of the hadronic MQ
systems (like those mentioned before) and the characteristics of hadronic
models that describe reliably microquasar jet emissions (Romero et al., 2016;
Smponias & Kosmas, 2011, 2014; Vieyro & Romero, 2017). The last few decades,
from the investigation of a great number of MQs and X-ray emitting binary
systems, stellar-mass black holes have been observed, with masses in the
region of our interest determined mainly from the dynamics and structure
properties of their companion stars (Remillard & McClintock, 2006; Charles &
Coe, 2006; Orosz, 2003).

The Galactic binary system SS433, located 5.5 kpc from the Earth in Aquila
constellation, displays two mildly relativistic jets (with bulk velocity
$\upsilon_{b}\approx 0.26$c) that are oppositely directed and precess in cones
(Romney et al., 1987). This system consists of a compact object (black hole)
and an A-type companion (donor) star in coupled orbit with a period $P\sim
13.1$ days (Fabrika, 2004). For the masses of the component stars of this
system, in this work we adopted the observations of the INTEGRAL
(Cherepashchuk et al., 2005), i.e. $M_{BH}=9M_{\odot}$ and
$M_{don}=30M_{\odot}$ for the black hole and the companion star, respectively.

Table 1: Parameters of the model for M33 X-7, Cygnus X-1 and SS433 Parameter | Symbol | SS433 | Cygnus X-1 | M33 X-7
---|---|---|---|---
Black Hole mass | $M_{BH}$ | $9.0M_{\odot}$ | $14.8M_{\odot}$ | $15.65M_{\odot}$
Distance from Earth | $d$ | 5.5 kpc | 1.86 kpc | 840-960 kpc
Donor Star mass | $M_{don}$ | $30M_{\odot}$ | $19.2M_{\odot}$ | $70M_{\odot}$
Donor star type | - | A-type | O-type | O-type
Orbital Period | $P$ | 13.1 days | 5.6 days | 3.45 days
Jet’s kinetic power | $L_{k}$ |  $10^{39}erg$ $s^{-1}$ |  $10^{38}erg$ $s^{-1}$ |  $10^{38}erg$ $s^{-1}$
Jet’s launching point | $z_{0}$ |  $1.3\times 10^{9}$ cm |  $10^{8}$ cm |  $1.9\times 10^{8}$ cm
Bulk velocity of jet particles | $\upsilon_{b}$ | 0.26c | 0.6c | 0.8c
Jet’s bulk Lorentz factor | $\Gamma_{b}$ | 1.04 | 1.25 | 1.66
Jet’s half-opening angle | $\xi$ | $0.6^{\circ}$ | $1.5^{\circ}$ | $7^{\circ}$
Jet’s viewing angle | $\theta$ | $78.05^{\circ}$ | $27.1^{\circ}$ | $74.6^{\circ}$

The second Galactic system addressed in this work, Cygnus X-1, is an X-ray
source in the constellation Cygnus which is located 1.86 kpc from the Earth
(Reid et al., 2011). The $19.2M_{\odot}$ O-type companion star (Orosz et al.,
2011) is in coupled orbit with the $14.8M_{\odot}$ compact object (stellar
mass black hole) (Orosz et al., 2011), with orbital period $P\sim 5.6$ days
(Gies et al., 2008).

As mentioned before, one of our main goals in this work is to calculate high
energy $\gamma$-ray and neutrino emissions from the recently discovered binary
system M33 X-7 (Long et al., 1981) which is located in the nearby galaxy
Messier 33, the only known up to now black hole that is in an eclipsing binary
system, with orbital period $P\sim 3.45$ days (Pietsch et al., 2006). The
distance from Earth of this system is between $\sim 840$ kpc (Orosz et al.,
2008) and $\sim 960$ kpc (Bonanos, Stanek et al., 2006). In addition, the
black hole mass is $M_{BH}=(15.65\pm 1.45)M_{\odot}$ which orbits its
$M_{don}=(70.0\pm 6.9)M_{\odot}$ O-type companion star.

In Table 1, we tabulate in addition some other important parameters employed
in this work for the systems under investigation, SS 433, Cygnus X-1 and M33
X-7.

###  2.1  Dynamics and energetic evolution of hadronic MQ jets

In the model considered in this work for the description of MQ jets, an
accretion disk in the equatorial region of the compact object is present, and
a fraction of the accreted material is expelled in two oppositely directed
jets (Falcke & Biermann, 1995; Smponias & Kosmas, 2011, 2014). We assume an
approximately conical mass outflow (jet) with a half-opening angle $\xi$ (for
the M33 X-7 MQ, we adopt the value $\xi=7^{\circ}$) and a radius given by
$r(z)=z\tan\xi$ (the coordinate z-axis coincides with the cone axis, jet’s
direction). The injection point (plane) of the jet is at a distance $z_{0}$
from the compact object. When $z=z_{0}$ the radius of the jet is given by
$r_{0}=z_{0}\tan\xi$ (see Table I).

The initial jet radius is $r_{0}=r(z_{0})\approx 5R_{sch}$, where
$R_{sch}=2GM_{BH}/c^{2}$ ($R_{sch}$ is the known Schwarzschild radius, the
critical black hole radius for which the escape speed of particles becomes
equal to the speed of light $c$), with $G$ denoting the Newton’s gravitational
constant. Then, we find, for example, that the injection point is at
$z_{0}=r_{0}/\tan\xi\simeq 1.9\times 10^{8}$ cm, for the M33 X-7. For the sake
of comparison, we mention that for Cygnus X-1, the injection point is at
$z_{0}\simeq 10^{8}$ cm, while for SS 433 it is $z_{0}\simeq 1.3\times 10^{9}$
cm.

It should be also noted that, for all systems studied in this work, the
extension of the jet is determined through a maximum value of z, denoted by
$z_{max}$, which is assumed to be equal to $z_{max}=5z_{0}$ (this point may be
assumed to be at the boundaries of the jet with the ambient region).

In discussing the energetic evolution and dynamics of MQ jets, the kinetic
energy density of the jet, $\rho_{k}(z)$, is related to its kinetic
luminosity, $L_{k}$, through the expression (Romero et al., 2016)

$\rho_{k}(z)=\frac{L_{k}}{\pi\upsilon_{b}{[r(z)]}^{2}}\,,$ (1)

where $\upsilon_{b}$ is the bulk velocity of the jet particles (mostly
protons). Within the context of the jet-accretion coupling hypothesis, only
around $10\%$ of the Eddington luminosity goes into the jet (Körding, Fender &
Migliari, 2006). Here we adopt $L_{k}=10^{38}erg\ s^{-1}$, for a black hole of
mass $\sim 15.7M_{\odot}$, i.e. for the M33 X-7 system.

Furthermore, assuming equipartition between the magnetic energy and the
kinetic energy in the jet, it implies that $\rho_{mag}=\rho_{k}$ (Bosch-Ramon,
Romero & Paredes, 2006), and hence the magnetic field that colimates the
astrophysical jet-plasma is given by (Reynoso & Romero, 2009)

$B(z)\equiv\sqrt{8\pi\rho_{mag}(z)}=\sqrt{8\pi\rho_{k}(z)}.$ (2)

In general, for the kinetic power in the jet, many authors consider that a
fraction is carried by primary protons and the rest by electrons, i.e.,
$L=L_{p}+L_{e}$. Then, the relation between the proton and electron power is
determined through a parameter $\alpha$ in such a way that $L_{p}=\alpha
L_{e}$ (Reynoso, Romero & Christiansen, 2008; Reynoso & Romero, 2009; Vieyro &
Romero, 2017).

In this article, however, we adopt the case of $\alpha=100$ which means that
we consider proton-dominated jet. Other quantities needed for the purposes of
this paper are listed in Table I.

###  2.2 p-p Collision Mechanism inside MQ jets

The collision of relativistic protons with the cold ones inside the jet (p-p
collision mechanism), produces high energy charged particles (pions
$\pi^{\pm}$, kaons $K^{\pm}$, muons $\mu^{\pm}$, etc.) and neutral particles
($\pi^{0}$, $K^{0}$, ${\tilde{K}}^{0}$, $\eta$ particles, etc.). Important
primary reactions of this type are

$\begin{split}p+p&\rightarrow p+p+a\pi^{0}+b(\pi^{+}+\pi^{-})\,,\\\
p+p&\rightarrow p+n+\pi^{+}+a\pi^{0}+b(\pi^{+}+\pi^{-})\,,\\\ p+p&\rightarrow
n+n+2\pi^{+}+a\pi^{0}+b(\pi^{+}+\pi^{-})\,,\end{split}$

where $a$ and $b$ denote the pion multiplicities (Romero et al., 2016), and
similarly for kaon production (Lipari, Lusignoli & Meloni, 2009).

Charged pions $\pi^{\pm}$ (and kaons $K^{\pm}$), afterwards, decay to charged
leptons (muons $\mu^{\pm}$, and electrons or positrons, $e^{\pm}$) as well as
neutrinos and anti-neutrinos as

$\begin{split}&\pi^{+}\rightarrow\mu^{+}+\nu_{\mu},\quad\quad\pi^{+}\rightarrow
e^{+}+\nu_{e}\\\
&\pi^{-}\rightarrow\mu^{-}+\bar{\nu_{\mu}},\quad\quad\pi^{-}\rightarrow
e^{-}+\bar{\nu_{e}}\,.\end{split}$ (3)

Furthermore, muons also decay giving neutrinos and electrons (or positrons) as

$\mu^{+}\rightarrow e^{+}+\nu_{e}+\bar{\nu_{\mu}},\quad\quad\mu^{-}\rightarrow
e^{-}+\bar{\nu_{e}}+\nu_{\mu}\,.$ (4)

We mention that, neutral pions ($\pi^{0}$), $\eta$-particles, etc. decay
producing $\gamma$-rays according to the reactions

$\pi^{0}\rightarrow\gamma+\gamma\,,\quad\quad\pi^{0}\rightarrow\gamma+e^{-}+e^{+}$
(5)
$\eta\rightarrow\gamma+\gamma\,,\quad\quad\eta\rightarrow\pi^{0}+e^{-}+e^{+}\,.$
(6)

In our present study, we consider neutrinos generated from both types of
charged-pion decays (for neutrinos coming from kaon decays the reader is
referred, e.g. to Ref. (Lipari, Lusignoli & Meloni, 2009) and also from both
types of charged muon decays (Papadopoulos, 2020; Papadopoulos, Papavasileiou
& Kosmas, 2020). Moreover, we consider $\gamma$-rays produced from the
$\pi^{0}$-decays of Eq. (5). Thus, high-energy neutrinos are inevitably
accompanied by pionic gamma-rays, a phenomenon well known as multi-messenger
emission from microquasar jets.

###  2.3  Accelerating and cooling rates of the jet processes

In general, the primary charged particles (p and $e^{-}$) gain energy during
moving within the magnetic field $B$ (Fermi acceleration). The acceleration
rate of the (initially cold) protons to an energy $E=E_{p}$, known as shock
acceleration, is defined as $t_{acc}^{-1}=E^{-1}dE/dt$ and is given by the
relation

$t_{acc}^{-1}\equiv\frac{1}{E}\frac{dE}{dt}\approx\eta\frac{ceB}{E_{p}}$ (7)

($\eta=0.1$ is the acceleration efficiency which means that only 10% of the
cold/thermal protons may acquire relativistic energies). This is equivalent to
existence of an efficient accelerator at the base of the jets where shocks are
rather relativistic (Begelman et al., 1980).

In the p-p collision mechanism, the cross section and the rate of inelastic
p-p scattering (relativistic protons scatter on cold ones),
$\sigma_{pp}^{inel}$ and $t_{pp}^{-1}$ respectively, play crucial role (the
in-elasticity coefficient is $K_{pp}\approx 1/2$). The corresponding cross
section, $\sigma_{pp}^{inel}$, as a function of the fast proton energy
$E_{p}$, reads (Kelner, Aharonian & Bugayov, 2006, 2009)

$\displaystyle\sigma_{pp}^{inel}(E_{p})=(0.25L^{2}+1.88L+34.3)$
$\displaystyle\left[1-\left(\frac{E_{th}}{E_{p}}\right)^{4}\right]^{2}$ (8)
$\displaystyle\times 10^{-27}cm^{2}$

where $L=ln(\frac{E_{p}}{1000\quad GeV})$, and $E_{th}$ denotes the minimum
energy of (thermal) protons, $E_{p}^{min}$, which is equal to
$E_{p}^{min}\equiv E_{th}=1.22$ GeV, see e.g. Appendix of Ref. (Smponias &
Kosmas, 2014).

The corresponding rate, $t_{pp}^{-1}$, of the inelastic p-p scattering is
given in Table 2 in terms of the cross section $\sigma_{pp}^{inel}$ and the
number density of cold particles (jet protons), $n_{p}(z)$, at a distance $z$
from the black hole which may be written as

$n_{p}(z)=\frac{(1-q_{rel})}{\Gamma_{b}m_{p}c^{2}}\rho_{k}(z)\,,$ (9)

($q_{rel}$=0.1 denotes the portion of the relativistic protons).

In Table 2, in addition to the p-p collision rate and accelerating rates of
protons, we tabulate the most significant cooling rates
(synchrotron,adiabatic, etc.) for protons, pions and muons as well.

Table 2: Accelerating, cooling and decay rates for particles moving inside hadronic MQ jets Rate Parameter | Rate Symbol | Basic rate definition
---|---|---
Proton accelerating rates | $t_{acc}^{-1}$ | $\eta\frac{ceB}{E_{p}}$
Proton-proton collision rate | $t_{pp}^{-1}$ | $n_{p}(z)\sigma_{pp}^{inel}(E_{p})K_{pp}$
Synchrotron radiation rate | $t_{syn}^{-1}$ | $\frac{4}{3}{\left(\frac{m_{e}}{m}\right)}^{3}\frac{\gamma\sigma_{T}B^{2}}{m_{e}c8\pi}$
Adiabatic expansion rate | $t_{ad}^{-1}$ | $\frac{2}{3}\frac{\upsilon_{b}}{z}$
Pion-proton collision | $t_{\pi p}^{-1}$ | $0.5\ c\ n_{p}(z)\sigma_{\pi p}^{inel}(E_{p})$
Proton escape rate | $t_{esc}^{-1}$ | $c/(z_{max}-z)$
Pion decay rate | $t_{\pi}^{-1}$ | $t_{esc,\pi}^{-1}+t_{dec,\pi}^{-1}$
Muon decay rate | $t_{\mu}^{-1}$ | $t_{esc,\mu}^{-1}+t_{dec,\mu}^{-1}$

In more detail, the well known synchrotron radiation is emitted from charged
particles of energy $E=\gamma mc^{2}$, with $\gamma$ being the particle’s
Lorentz factor and $m$ its mass, moving inside the magnetized plasma. Thus,
the rate of the synchrotron radiation, $t_{syn}^{-1}$, is very important in MQ
systems with strong magnetic fields (see Table 2).

Furthermore, due to the adiabatic expansion of the jet plasma, all co-moving
particles inside it loose energy with a rate $t_{ad}^{-1}$ (Bosch-Ramon,
Romero & Paredes, 2006). This rate, known as adiabatic cooling rate, depends
on the particle velocities $\upsilon_{b}\tan\xi$.

In addition, $t_{\pi p}^{-1}$ gives the rate of pion-proton ($\pi-p$)
inelastic scattering inside the jet. The corresponding $\sigma_{\pi p}^{inel}$
cross section is related to the inelastic p-p scattering cross section with
the relationship $\sigma_{\pi p}(E)=(2/3)\sigma_{pp}^{inel}(E)$, where the
factor $2/3$ comes out of the fact that protons are made of three valence
quarks, while the pions of only two quarks (Gaisser, 1990).

It is worth mentioning that some processes taking place inside the jet, lead
to particles’ knock out, either by escaping from the jet (and/or from the
energy range of interest) or by their decay processes (Ginzburg & Syrovatskii,
1964). Thus, the corresponding rates are: (i) the escaping rate,
$t^{-1}_{esc}$, given by

$t_{esc}^{-1}\approx\frac{c}{z_{max}-z}\,,$ (10)

(${z_{max}-z_{0}}$ represents the length of the acceleration zone), and (ii)
the decay rate of the particle in question, $t^{-1}_{dec}$, which is known
from direct measurements of particle’s life time.

The rates of the knock out (or catastrophic) processes inside the jet for
pions, $t_{\pi}^{-1}$, and muons, $t_{\mu}^{-1}$, include contributions from
the decay process and escaping process as

$t_{\pi}^{-1}(E,z)=t_{esc,\pi}^{-1}(z)+t_{dec,\pi}^{-1}(E)\,,t_{\mu}^{-1}(E,z)=t_{esc,\mu}^{-1}(z)+t_{dec,\mu}^{-1}(E)\,,$
(11)

where $t_{dec,\pi}^{-1}={[(2.6\times 10^{-8})\gamma_{\pi}]}^{-1}\ s^{-1}$, for
pions, and $t_{dec,\mu}^{-1}={[(2.2\times 10^{-6})\gamma_{\mu}]}^{-1}\
s^{-1}$, for muons. We note that Eq. (11) holds also for protons due to the
fact that $t^{-1}_{dec,p}=0$ (or equivalently the proton life time is
infinite).

The rates discussed above enter the set of basic coupled transport equations
(system of kinematic equations) that describe the particle distributions as we
will discuss in Section 3.

####  2.3.1  Energy variation of the various rates of jet particles

In Figs. 1, 2 and 3, the variation of the aforementioned main rates, versus
the particle energy E, for protons, pions and muons, respectively, and for the
SS 433 and M33 X-7 MQ systems, are illustrated. In other words, in these
figures the region of dominance of the accelerating, cooling, etc. rates
throughout the energy range of interest, $E_{th}\equiv E^{min}\leq E\leq
E^{max}$, are demonstrated.

Figure 1:  Cooling rates for protons in the jets of M33 X-7 (top) and SS433
(bottom) at the base of the jets $z_{0}$. The plots of cooling rates show the
synchrotron emission (solid lines), the adiabatic cooling (dotted lines), and
the inelastic p-p collision rate (dashed lines).

More specifically, Fig. 1 indicate the cooling rates at the base of the jet
for protons, Fig. 2 those for pions and Fig. 3 the corresponding ones for
muons, in the extragalactic binary system M33 X-7 (top of these figures). For
comparison, the corresponding results for the Galactic binary system SS 433
(bottom of these figures) are also shown.

The plots in Fig. 2 for pions, show respectively the variation versus the
energy $E$ of the rate for synchrotron emission (solid lines), for the pion-
proton inelastic collision (dashed lines), for the adiabatic cooling (dotted
lines), and for the decay of pions (dot-dashed lines). They refer to the M33
X-7 (top) and the SS 433 (bottom) micro-quasar systems. Also, Fig. 3
illustrates the cooling rates of muons for the synchrotron emission (solid
lines), the adiabatic expansion (dotted lines), and for the decay of muons
(dot-dashed lines). They are similar to those for pions (Fig. 2) except the
corresponding muon-proton and muon-pion scattering channels which are ignored
in this work.

Figure 2: Cooling rates for pions in the jets of M33 X-7 (top) and SS433
(bottom) at the base of the jets

From Figs. 2 and 3 it becomes obvious that, the particle synchrotron losses
dominate the high energy region. In the case of protons (Fig. 1), however, due
to the large proton mass, synchrotron losses are not dominant up to very high
energies ($\sim 10^{7}$ GeV) for the MQ systems studied. On the other hand,
for pions and muons, the decay losses dominate for lower energies, due to
their short life time (see Table 2).

The main difference between these two MQ systems (SS 433 and M33 X-7) concerns
the synchrotron cooling rates. In general, according to Eqs. (1) and (2), for
systems with wide half-opening angle (e.g. M33 X-7, with $\xi=7^{\circ}$), the
magnetic energy density is lower than that of systems with small $\xi$ (e.g.
the SS 433 with $\xi=0.6^{\circ}$). Hence, the magnetic field for M33 X-7 is
lower which becomes obvious by comparing Eq. (2) for the two systems. The
latter conclusion justifies the lower synchrotron loss rate for the M33 X-7
system.

Figure 3: Cooling rates for muons in the jets of M33 X-7 (top) and SS433
(bottom) at the base of the jets

##  3  Formulation of the key-role transport equations inside jet plasmas

In this section, we make a first step towards improving the mathematical
formulation describing the creation, scattering, decay, diffusion and emission
of particles (including neutrinos) and electromagnetic radiation in
astrophysical outflows (jets). The basic key-role tool of this formalism is
the general transport equation which satisfies each of the particles
considered (Vieyro & Romero, 2017).

We should note that, such a formalism may be reliably applicable for multi-
species (multi-particle) emissions from MQs and X-ray binaries, as well as
from central regions of galaxies like the AGN systems, in which the compact
object is a massive or supermassive black hole, see e.g. (Romero et al., 2016;
Ginzburg & Syrovatskii, 1964). The geometry of the latter sources is assumed
to be rather spherical, while in the first class of systems that are studied
in this work, the sources are assumed conical and the particles (neutrinos) or
radiation are emitted to the forward or backward direction.

From a mathematical point of view, the transport equation which satisfies any
kind of particles moving into the (dark) jets of MQs, is an integro-
differential equation, as is explained below.

###  3.1 The transport equation for particles moving inside MQ jets

The general transport equation describes the concentration (distribution) of
particle of j-kind, $N_{j}(E,{\bf r},t)$, where
$j=p,e^{\pm},\pi^{\pm},\mu^{\pm}$, etc., as a function of the time $t$, the
particle’s energy $E$, and the position ${\bf r}$ inside the (conical) jet. In
essence, this is a phenomenological macroscopic equation of particle (or
radiation) transport describing astrophysical outflows, and it is written as

$\frac{\partial N_{j}}{\partial t}-\nabla\cdot\left(D_{j}\nabla
N_{j}\right)+\frac{\partial\left(b_{j}N_{j}\right)}{\partial
E}-\frac{1}{2}\frac{\partial^{2}\left(d_{j}N_{j}\right)}{\partial
E^{2}}=Q_{j}(E,{\bf r},t)-p_{j}N_{j}+\sum_{k}\int
P_{j}^{k}(E^{\prime},E)N_{j}(E,{\bf r},t)dE\,,\quad
j=p,\pi^{\pm},\mu^{\pm}\,...$ (12)

where the parameters $D_{j}$, $b_{j}$, $d_{j}$, $p_{j}$ and $P_{j}^{k}$ may
depend upon the space and time coordinates and also on the energy $E$. The
latter equation coincides with the continuity equation for particles of
type-$j$, with j=1,2,3 and $(1,2,3)\equiv$(p, $\pi$, $\mu$)].

The term $Q_{j}(E,{\bf r},t)$ in the r.h.s. of Eq. (12) is equal to the
intensity of the source producing the particles-$j$, which is also known as
the injection function of particles-j. This means that $Q_{j}(E,{\bf
r},t)dEd^{3}{\bf r}dt$ represents the number of particles kind-$j$ provided by
the sources in a volume element $d^{3}{\bf r}$, in the energy range between
$E$ and $E+dE$ during the time $dt$. In the case when the $j$-type particles
are products of a chain reaction (as it holds in our present work assuming the
p-p reaction chain described in Sect. 2.2), the function $Q_{j}(E,{\bf r},t)$
couples the $j$-reaction with its parent reaction, i.e. the equation of
particles kind-(j-1).

Further, the term proportional to $p_{j}$ exists in cases when catastrophic
(or knock out) processes take place that cause catastrophic energy losses.
Then, this gives the probability per unit time for the losses in question to
occur. Thus, the coefficient $p_{j}$ is written as

$p_{j}=\frac{1}{T_{j}}\equiv t^{-1}_{j}\,,$ (13)

with $T_{j}$ is the mean life time of the particles of kind-$j$. Also
$p_{j}N_{j}$ denotes the number of particles "knocked out" per unit time.

The coefficients $b_{j}$ in the third term (r.h.s.) are equal to the mean
energy increment of the particle-$j$ per unit time i.e.

$b_{j}=\frac{dE}{dt}\,.$ (14)

For the special case when $b_{j}=b_{j}(E)$, this coefficient is related to
energy losses of various cooling processes. Also, the coefficients
$P^{k}_{j}(E)$, with $k>j$, in the latter summation of the transport equation,
are related with the existence of fragmentation of the primary (or secondary)
particles which in this work are assumed as non-existing.

The parameter $D_{j}$ is known as the diffusion coefficient which, in general,
is a functions of the coordinates $\bf r$ and time $t$ if the concentration of
the jet plasma and its macroscopic motion are inhomogeneous in the volume of
the jet. By assuming that $D_{j}=0$ we may reliably describe the regular
motion of particles j-kind along the lines of force of the magnetic field
${\bf B}$.

For simplicity, in this work we make the realistic assumptions that, the above
coefficients depend only on the particle’s energy $E$, i.e. we assume that
$D_{j}(E)$, $b_{j}(E)$ and $p_{j}(E)$. Then, the system of general transport
equations (12) reads

$\frac{\partial}{\partial t}\begin{pmatrix}N_{p}(E,z,t)\\\ N_{\pi}(E,z,t)\\\
N_{\mu}(E,z,t)\end{pmatrix}-\nabla^{2}\begin{pmatrix}D_{p}(E)N_{p}(E,z,t)\\\
D_{\pi}(E)N_{\pi}(E,z,t)\\\
D_{\mu}(E)N_{\mu}(E,z,t)\end{pmatrix}+\frac{\partial}{\partial
E}\begin{pmatrix}b_{p}(E)N_{p}(E,z,t)\\\ b_{\pi}(E)N_{\pi}(E,z,t)\\\
b_{\mu}(E)N_{\mu}(E,z,t)\end{pmatrix}-\frac{1}{2}\frac{\partial^{2}}{\partial
E^{2}}\begin{pmatrix}d_{p}(E)N_{p}(E,z,t)\\\ d_{\pi}(E)N_{\pi}(E,z,t)\\\
d_{\mu}(E)N_{\mu}(E,z,t)\end{pmatrix}+\begin{pmatrix}t_{esc}^{-1}N_{p}(E,z,t)\\\
t_{\pi}^{-1}N_{\pi}(E,z,t)\\\
t_{\mu}^{-1}N_{\mu}(E,z,t)\end{pmatrix}=\begin{pmatrix}Q_{p}(E,z,t)\\\
Q_{\pi}(E,z,t)\\\ Q_{\mu}(E,z,t)\end{pmatrix}$ (15)

Note that, in the latter system of coupled differential equations we write
down only three particles (protons, pions and muons), without considering the
complete reaction family tree, since we haven’t distinguished particles with
different charge as $\pi^{\pm}$ and $\mu^{\pm}$, and we haven’t considered the
possibility of left-right symmetry, $\mu_{L,R}$ of muons as we have done below
in Sect. 4. This means that, by considering all these particles, the system of
Eq. (15) will have seven lines.

Usually, we ignore the term involving second derivative with respect to $E$
means that the acceleration term is negligible.

The solution of the general transport equation for particles of one kind (when
the last term of this equation can be omitted) is described in detail in Ref.
(Kosmas & Smponias, 2018). In the latter work, however, neutrinos produced
only from $\pi^{\pm}$ have been considered. In the present paper we proceed
further and include neutrino emissions also from the muon decays, so we need
to solve the system of Eq. (15) by including the third line too.

The complete form of the system of Eq. (15) is treated with the method of Ref.
(Tsoulos, Kosmas & Stavrou, 2017) in order to find the exact solutions. In
Ref. (Papadopoulos, Papavasileiou & Kosmas, 2020), equation (15) is treated
semi-analytically. In the present work, we restrict ourselves to simplified
forms resulting by neglecting various phenomena (processes) taking place
inside the astrophysical outflows (jet plasma), even though some of the
assumed omissions may be considered as rather crud approximations.

Below we discuss the cases of Eq. (15) satisfying the conditions of the one-
zone approximation (Khangulyan et al., 2007; Kosmas & Smponias, 2018), i.e.
the cases when the particle distributions are independent of time (steady
state approximation).

####  3.1.1  Transport equation assuming absence of energy-losses

At first, in the calculations of our present work, we start by writing the
simplest solutions of the system of transport equations obtained by assuming
that all energy losses are absent i.e. $b_{j}=0$. Then, by denoting $N_{p,0}$,
$N_{\pi,0}$ and $N_{\mu,0}$ the corresponding energy distributions for
protons, pions and muons, respectively, the system of transport equations
takes the trivial form

$\begin{pmatrix}t_{p}^{-1}(E)N_{p,0}(E,z)\\\ t_{\pi}^{-1}(E)N_{\pi,0}(E,z)\\\
t_{\mu}^{-1}(E)N_{\mu,0}(E,z)\end{pmatrix}=\begin{pmatrix}Q_{p}(E,z)\\\
Q_{\pi}(E,z)\\\ Q_{\mu}(E,z)\end{pmatrix}$ (16)

The calculated distributions $N_{j,0}$, with j=p, $\pi^{\pm}$, $\mu^{\pm}$ of
the latter equations are discussed in the next Section.

####  3.1.2 Steady state transport equations with particle losses when moving
inside the jet plasma

Under the conditions of the steady state approximation, and assuming that the
various energy losses are absent, the system of transport equations Eq. (15)
takes the form

$\frac{\partial}{\partial E}\begin{pmatrix}b_{p}(E)N_{p}(E,z)\\\
b_{\pi}(E)N_{\pi}(E,z)\\\
b_{\mu}(E)N_{\mu}(E,z)\end{pmatrix}+\begin{pmatrix}t_{esc}^{-1}N_{p}(E,z)\\\
t_{\pi}^{-1}N_{\pi}(E,z)\\\
t_{\mu}^{-1}N_{\mu}(E,z)\end{pmatrix}=\begin{pmatrix}Q_{p}(E,z)\\\
Q_{\pi}(E,z)\\\ Q_{\mu}(E,z)\end{pmatrix}$ (17)

In order to calculate the time independent neutrino and gamma-ray
emissivities, we need, first, to calculate the distributions of protons, pions
and muons, $N_{p}$, $N_{\pi}$, $N_{\mu}$, respectively, from Eq. (17). For
these computations we used a code written in the C programming language,
mainly following the assumptions of Refs. (Reynoso, Romero & Christiansen,
2008; Reynoso & Romero, 2009; Vieyro & Romero, 2017).

In the latter case, the second term of the l.h.s. of the transport equation we
replace the escape rate, $t_{esc}^{-1}$ of Eq. (10) and
$b(E)=-Et_{loss}^{-1}(E)\,$. The latter quantity the energy loss rate for each
particle which is given by

$\displaystyle b_{p}(E)$
$\displaystyle=-Et_{loss,p}^{-1}(E)=-E(t_{syn}^{-1}+t_{ad}^{-1}+t_{pp}^{-1}),$
(18) $\displaystyle b_{\pi}(E)$
$\displaystyle=-Et_{loss,\pi}^{-1}(E)=-E(t_{syn}^{-1}+t_{ad}^{-1}+t_{\pi
p}^{-1}),$ (19) $\displaystyle b_{\mu}(E)$
$\displaystyle=-Et_{loss,\mu}^{-1}(E)=-E(t_{syn}^{-1}+t_{ad}^{-1})$ (20)

It is worth noting that, in the latter equations some additional terms may
appear in cases when the rates of some other processes, ignored in our present
work (like e.g. the Inverse Compton Scattering, the pion-muon scattering,
etc.), may be considered important for the description of other cosmic
structures.

##  4 Injection functions and energy distributions of particles inside the jet

In this section, the source functions $Q_{j}(E,{\bf r},t)$ entering the r.h.s
of the system of coupled integro-differential equations of transport type, are
discussed. Towards this aim, initially in the system of Eqs. (15), we insert
phenomenological expressions obtained from Ref. (Lipari, Lusignoli & Meloni,
2009). A mathematical semi-analytic way of solving this system is proposed in
Ref. (Kosmas, 2020). Also, following Ref. (Tsoulos, Kosmas & Stavrou, 2017), a
method is derived for the numerical solution of the differential system of
transport equations.

###  4.1 Proton injection function and energy distribution

####  4.1.1 Relativistic proton injection function

In the assumed p-p mechanism, the density of fast (relativistic) protons
injected is more important, while the corresponding density of slow (thermal)
protons is dynamically significant.

A usual injection function for the relativistic protons, coming out of the
acceleration mechanism, has been taken to be a power-law with exponent equal
to two, i.e. a function of the energy of the form (Ghisellini, Maraschi &
Treves, 1985)

$\displaystyle
Q_{p}(E,z)=Q_{0}\left(\frac{z_{0}}{z}\right)^{3}\frac{1}{E^{2}}\,,$ (21)

where $Q_{0}$ is a normalization constant obtained by specifying the power in
the relativistic protons (Reynoso, Romero & Christiansen, 2008), see Appendix.
The above form is valid for the jet’s frame of reference (the corresponding
expression in observer’s system is shown in Appendix).

It should be noted that, Eq. (21) enters (through integration) the
determination of the injection functions of the secondary particles
($\pi^{\pm}$, $\mu^{\pm}$, see below) and this justifies the term systems of
coupled integro-differential equations used in Sect. 3. The latter injection
functions are also dependent on the rates and cross sections of the reactions
that preceded in the chain processes following the inelastic p-p collision.

####  4.1.2 Relativistic Proton energy distribution in the jet

Under the assumptions discussed before, the system of transport equations (17)
can be easily solved and the obtained proton distribution $N_{p}(E,z)$ is
written as (Reynoso, Romero & Christiansen, 2008; Romero & Vila, 2008)

$N_{p}(E,z)=\int_{E_{p}}^{E_{p}^{max}}\mid
b_{p}(E)\mid^{-1}Q_{p}(E^{\prime},z)e^{-\tau_{p}(E,E^{\prime})}dE^{\prime}$
(22)

where

$\tau_{p}(E,E^{\prime})=\int_{E}^{E^{\prime}}\frac{t_{esc,p}^{-1}dE^{\prime\prime}}{\mid
b_{p}(E^{\prime\prime})\mid}\,.$ (23)

The quantity $Q_{p}(E,z)$ corresponds to the relativistic injection function
of protons at the observer’s frame (see Appendix). The minimum energy of
protons is $E_{p}^{min}=E_{th}=1.22$ GeV, while the maximum energy is assumed
to be $E_{p}^{max}=10^{6}-10^{7}$ GeV.

Figure 4: Proton energy distribution for M33 X-7 and SS 433 microquasar
systems.

Figure 4 shows the proton energydistribution $N_{p}(E)$ at the base of the jet
for M33 X-7 (solid line) and SS 433 (dashed line). As can be seen, the number
of protons appears reduced (even by two orders of magnitude) for the wide
half-opening angle system of M33 X-7 ($\xi=7^{\circ}$) as compared to the
narrow of SS 433 ($\xi=0.6^{\circ}$). This is because, the jet radius $r(z)$
is bigger and hence the magnetic field $B$ is smaller in the case of M33 X-7
than that of SS433.

###  4.2 Pion injection functions and energy distributions inside the jet

####  4.2.1 Pion injection functions

For pions produced through the inelastic p-p scattering, as injection function
we adopt the one given by (Kelner, Aharonian & Bugayov, 2006) as

$\displaystyle
Q_{\pi}(E,z)=n(z)c\int_{\varepsilon}^{1}\sigma_{pp}^{inel}\left(\frac{E}{x}\right)N_{p}\left(\frac{E}{x},z\right)F_{\pi}\left(x,\frac{E}{x}\right)\frac{dx}{x}$
(24)

where $\varepsilon=E/E_{p}^{max}$, $x=E/E_{p}$ and $F_{\pi}$ denotes the
distribution of pions produced per $p-p$ collision (see Appendix).

####  4.2.2 Pion energy distributions

The steady state energy distribution $N_{\pi}(E,z)$ for pions (created through
the scattering of hot protons off the cold ones), results through a similar
way to that of Eqs. (22) and (23). With the replacement $p\rightarrow\pi$ and
$t_{esc}^{-1}\rightarrow t_{\pi}^{-1}(E,z)$, respectively, on the latter
equations the corresponding solution for pion energy distribution is, then,
written as

$\displaystyle N_{\pi}(E,z)$ $\displaystyle=\int_{E}^{E_{max}}\mid
b_{\pi}(E)\mid^{-1}Q_{\pi}(E^{\prime},z)e^{-\tau_{\pi}}dE^{\prime}$ (25)

where the rate $t_{\pi}^{-1}$ includes contributions from the decay and escape
rates [see Eq. (11)].

###  4.3 Muon injection functions and energy distributions inside the jet

####  4.3.1 Muon injection functions

As mentioned before, in this work in addition to neutrinos coming from the
decay of charged pions $\pi^{\pm}$, we evaluate neutrino emissivities
originating from the decay of the secondary charged muons $\mu^{\pm}$. For the
latter emissivity, we follow Ref. (Lipari, Lusignoli & Meloni, 2009), in order
to take into account properly the muon energy loss. This means that, it is
necessary to consider the production of both left handed and right handed
muons, $\mu_{L}^{-}$ and $\mu_{R}^{+}$, separately, because $\mu_{L}^{-}$ and
$\mu_{R}^{+}$ have different decay spectra.

Thus, the injection functions of the left handed and right handed muons are
(Lipari, Lusignoli & Meloni, 2009):

$\displaystyle
Q_{\mu_{L}^{-},\mu_{R}^{+}}(E_{\mu},z)=\int_{E_{\mu}}^{E^{max}}$
$\displaystyle t_{dec,\pi}^{-1}(E_{\pi})N_{\pi}(E_{\pi},z)\Theta(x-r_{\pi})$
(26)
$\displaystyle\times\frac{r_{\pi}(1-x)}{E_{\pi}x{(1-r_{\pi})}^{2}}dE_{\pi}$

and

$\displaystyle
Q_{\mu_{R}^{-},\mu_{L}^{+}}(E_{\mu},z)=\int_{E_{\mu}}^{E^{max}}$
$\displaystyle t_{dec,\pi}^{-1}(E_{\pi})N_{\pi}(E_{\pi},z)\Theta(x-r_{\pi})$
(27)
$\displaystyle\times\frac{(x-r_{\pi})}{E_{\pi}x{(1-r_{\pi})}^{2}}dE_{\pi}$

with $x=E_{\mu}/E_{\pi}$ and $r_{\pi}={(m_{\mu}/m_{\pi})}^{2}$.

####  4.3.2 Muon energy distributions

The corresponding muon energy distributions, coming out of the transport
equation, result from Eq. (25) via the replacement $t_{esc}^{-1}\rightarrow
t_{\mu}^{-1}(E,z)$, for muons and are then written as

$\displaystyle N_{\mu_{i}}(E,z)$ $\displaystyle=\int_{E}^{E_{max}}\mid
b_{\mu}(E)\mid^{-1}Q_{\mu_{i}}(E^{\prime},z)e^{-\tau_{\mu}}dE^{\prime}$ (28)

where the rate $t_{\mu}^{-1}$ for muons is given in Eq. (11).

Figure 5: Pion (Top) and muon (bottom) distributions for M33 X-7. In both case
of particles-j, we compare the distribution which considers the important
energy losses $N_{j}(E,z=z_{0})$ and with that which neglects energy losses
$N_{j,0}(E,z=z_{0})$ for $j=\pi,\ \mu$. Both types of particle distributions
refer to the base of the jet, i.e. at $z=z_{0}$

In Fig. 5, the pion (top) and muon (bottom) energy distributions for M33 X-7
are illustrated. In both cases of particles, the distribution
$N_{j}(E,z=z_{0})$, with $j=\pi$ or $\mu$, which considers the most important
energy losses, is compared with that which neglects energy losses,
$N_{j,0}(E,z=z_{0})$. The solid lines correspond to distributions obtained by
considering energy losses and the dashed lines correspond to those where the
energy losses obtained have been neglected.

##  5 Results for neutrino and $\gamma$-ray emissivities

In this Section, we present and discuss simulated emissivities of high energy
neutrinos and $\gamma$-rays produced in the extragalactic microquasar M33 X-7
system which is located in the Messier 33 galaxy, at a distance $\sim 840-960$
kpc from the Earth (Pietsch et al., 2006). The derived algorithms for this
purpose have been tested on the well-studied Galactic micro-quasars SS 433 and
Cyg X-1 system as well.

Because in our previous calculations on neutrino production from MQs (Smponias
& Kosmas, 2015, 2017; Kosmas & Smponias, 2018), we neglected (due to the
complexity and long time consuming) emissivity (intensity) of neutrinos
originated from the secondary muons ($\mu^{\pm}$) produced from the charged
pion ($\pi^{\pm}$) decays, in this section, we present contributions also
originating from the $\mu^{\pm}$ channel [see Eq. (4)].

Figure 6: Neutrino (left) and gamma-ray (right) intensity for M33 X-7 and
SS433

By using the concentrations for protons $N_{p}$, pions $N_{\pi}$ and muons
$N_{\mu}$, the neutrino and $\gamma$-ray intensities are subsequently
calculated as described below.

###  5.1  Neutrino emission from the p-p reaction family tree

After the above discussion, the total emissivity $Q_{\nu}(E,z)$ produced from
a MQ is the sum of contributions from the two sources: (i) the first comes
from the direct $\pi^{\pm}$ decay (prompt neutrino production), and (ii) the
second comes from the $\mu^{\pm}$ decay (delayed neutrino production). Thus,

$Q_{\nu}(E,z)=Q_{\pi\rightarrow\nu}(E,z)+Q_{\mu\rightarrow\nu}(E,z)$ (29)

For pion decays the injection function is given by

$\displaystyle
Q_{\pi\rightarrow\nu}(E,z)=\int_{E}^{E_{max}}t_{dec,\pi}^{-1}(E_{\pi})N_{\pi}(E_{\pi},z)\frac{\Theta(1-r_{\pi}-x)}{E_{\pi}(1-r_{\pi})}dE_{\pi}$
(30)

with $x=E/E_{\pi}$, while for the four types of muon decays the injection
function reads

$\displaystyle
Q_{\mu\rightarrow\nu}(E,z)=\sum_{i=1}^{4}\int_{E}^{E_{max}}t_{dec,\mu}^{-1}(E_{\mu})N_{\mu_{i}}(E_{\mu},z)$
$\displaystyle\times\left[\frac{5}{3}-3x^{2}+\frac{4}{3}x^{3}+\left(3x^{2}-\frac{1}{3}-\frac{8}{3}x^{3}\right)h_{i}\right]\frac{dE_{\mu}}{E_{\mu}}\,.$
(31)

(with $x=E/E_{\mu}$).

In the last expression, the symbols $\mu_{i},i=1,2,3,4$ correspond to
$\mu_{\\{1,2\\}}=\mu_{L}^{\\{-,+\\}}$, $\mu_{\\{3,4\\}}=\mu_{R}^{\\{-,+\\}}$
(Lipari, Lusignoli & Meloni, 2009), while $h_{\\{1,2\\}}=-h_{\\{3,4\\}}=-1$.
Then, the calculation of each of the latter two integrals of Eqs. (30) and
(31) provides separately the partial emissivity of neutrinos for the prompted
and the delayed neutrino source, respectively. Needless to note that Earth and
space telescope are not able to discriminate the two source as it happens with
laboratory neutrino sources.

Subsequently, we easily obtain the neutrino intensity (in units
$GeV^{-1}s^{-1}$) by the spatial integration

$I_{\nu}(E)=\int_{V}Q_{\nu}(E,z)d^{3}r=\pi\left(\tan\xi\right)^{2}\int_{z_{0}}^{z_{max}}Q_{\nu}(E,z)z^{2}dz\,.$
(32)

Figure 6 (left) shows the neutrino intensity produced at the base of the jet
from direct decays of secondary pions and muons coming from p-p collisions in
the jets of M33 X-7 and SS 433, respectively. As can be seen, the number of
produced neutrinos reduces significantly for energies $E>10^{3}$ GeV,
following the behavior of pion and muon distribution. We can also see that the
neutrinos produced in M33 X-7 are less than those produced in SS 433, due to
the wider-half opening angle $\xi$ as we explained before.

It should be noted that, to perform the above calculations we have updated and
improved our codes to reduce the time consuming to a reasonable level, so as
the integral involved in Eq. (31) to be easily obtained, since in going the
step from Eq. (30) to Eq. (31) the time consuming increases rapidly.

###  5.2  Gamma-ray emission from the p-p reaction chain

As we have discussed in Sect. 2, the p-p collisions inside the jets, produce
secondary neutral particles (pions, eta particle, etc.) that decay to give
$\gamma$-rays. The dominant of these channels goes through the reaction

$\pi^{0}\rightarrow\gamma+\gamma\,.$ (33)

For $E_{\gamma}\geq 100$ $GeV$, we consider the $\gamma$-ray emissivity at a
distance $z$ along the jets (in units $GeV^{-1}s^{-1}$) to be given by
(Reynoso, Romero & Christiansen, 2008) as

$Q_{\gamma}=c\int_{E_{\gamma}/E_{p}^{max}}^{1}\sigma_{pp}^{inel}\left(\frac{E_{\gamma}}{x}\right)N_{p}\left(\frac{E_{\gamma}}{x},z\right)F_{\gamma}\left(x,\frac{E_{\gamma}}{x}\right)\frac{dx}{x}\,.$
(34)

where $F_{\gamma}$ is the spectrum of the produced $\gamma$-rays (Kelner,
Aharonian & Bugayov, 2006, 2009) with energy $x=E_{\gamma}/E_{p}$ for a
primary proton energy $E_{p}$. The function $F_{\gamma}$ is given in Appendix.

Subsequently, the corresponding spectral intensity of $\gamma$-rays, can be
obtained from the following spatial integration over the jet’s volume $V$ as

$I_{\gamma}(E_{\gamma})=\int_{V}Q_{\gamma}(E_{\gamma},z)d^{3}r=\pi\left(\tan\xi\right)^{2}\int_{z_{0}}^{z_{max}}Q_{\gamma}(E_{\gamma},z)z^{2}dz,$
(35)

Figure 6 (right) shows the $\gamma$-ray intensity for energies
$E_{\gamma}>100$ GeV. As can be seen, the produced $\gamma-ray$ intensity
reduces rather steadily following the behavior of the proton distribution in
both cases of MQ systems. Again for the M33 X-7 it is lower mainly due to the
wider half-opening angle (i.e. weaker magnetic field) but also due to other
parameters.

In our ongoing calculations (Kosmas, 2020), among others we include results
obtained by using the 3-D PLUTO astrophysical hydrocode where the jets of the
studied systems are approximated as purely relativistic magnetohydrodynamical
(RMHD) flows with the magnetic field (assumed tangled) being dragged with the
flow (Smponias & Kosmas, 2014).

##  6 Summary and Conclusion

In this work we address neutrino and $\gamma$-ray emissions from microquasars
and X-ray binary stars (XRBs) that consist of a stellar mass black hole
(compact object) and a main sequence donor star. We assume that these
emissions originate from decay and scattering processes of the secondary
particles produced through the $p-p$ scattering mechanism, i.e. the inelastic
collision of relativistic protons of the jet with the thermal ones.

Such high energy neutrinos and $\gamma$-rays are detectable by operating
terrestrial and space telescopes. Among the operating detectors are the under
ice IceCube (at the South Pole), the ANTARES and KM3Net (under the
Mediterranean sea), and many others.

In the near future, more accurate measurements of gamma-ray and neutrino
fluxes by the coming observatories, such as the CTA (Acharya et al.2018),
which is expected to shed more light on the nature and the emission sources.
From the perspective of neutrino detection, the addition of more years of data
with continuous operation of IceCube will improve the sensitivity of the
search for Galactic sources of cosmic neutrinos. Furthermore, the next
generation detection instrument, IceCube-Gen2, a substantial expansion of
IceCube, will be 10 times larger. This next-generation neutrino observatory
with five times the effective area of IceCube is expected to improve the
neutrino source search sensitivity by the same order (Aartsen et al., 2019).
With higher neutrino statistics, identifying Galactic sources will become more
promising.

On the other hand, theoretically, by modelling the solution of the system of
coupled transport equations, we were able to perform detailed calculations for
various processes taking place inside the jets of Galactic and extra galactic
system M33 X-7 assuming hadronic content in their jets. These cooling rates
enter the proton, pion and muon energy distributions through which one obtains
neutrino and $\gamma$-ray intensities.

Our near future calculations include, among others, multidimentional
simulations using the PLUTO astrophysical hydrocode which may provide both
radiative and dynamical description of the astrophysical outflows within
multi-scale and multi-messesger astrophysics.

## Acknowledgements

This research is co-financed (O.T.K.) by Greece and the European Union
(European Social Fund-ESF) through the Operational Programme “Human Resources
Development, Education and Lifelong Learning 2014- 2020" in the context of the
project (MIS5047635). D.A.P. wishes to thank Prof. T.S. Kosmas for fruitful
discussions during his stay in the Dept. of Physics, University of Ioannina.

## Data Availability

There is no data used in this paper.

## References

  * Aartsen et al. (2015) Aartsen M. G. et al., 2015 ApJ, 805, L5
  * Aartsen et al. (2016) Aartsen M. G. et al., 2016, (IceCube) Phys.Rev.Lett., 117 no.7, 071801.
  * Aartsen et al. (2018) Aartsen M. G. et al., 2018, (IceCube) Science, 361, 147
  * Aartsen et al. (2019) Aartsen M. G. et al., 2019, (IceCube) arXiv:1911.02561v1 [astro-ph.HE]
  * Aartsen et al. (2020) Aartsen M. G. et al., 2020, (IceCube Collaborations), Phys. Rev. Lett., 124, 051103, arXiv:1910.08488 [astro-ph.HE].
  * Actis et al. (2011) Actis M., Agnetta G., Aharonian F. et al., 2011, Exp. Astron., 32, 193
  * Adrian-Martinez et al. (2016) Adrian-Martinez S. et al., 2016, (KM3Net Collaboration) J. Phys. G: Nucl. Part. Phys., 43, 084001, e-Print: arXiv: [astro-ph.IM] 1601.07459
  * Albert et al. (2020) Albert A. et al., 2020, (ANTARES and IceCube Collaborations), ApJ, 892, 9, arXiv:2001.04412 [astro-ph.HE]
  * Begelman et al. (1980) Begelman M. C., Hatchett S. P.,McKee C. F., Sarazin C. L., Arons J., 1980, ApJ, 238, 722
  * Bonanos, Stanek et al. (2006) Bonanos A. Z., Stanek K. Z. et al., 2006, ApJ, 652, 313
  * Bosch-Ramon, Romero & Paredes (2006) Bosch-Ramon V., Romero G. E., Paredes J. M., 2006, A&A, 447, 263
  * Charles & Coe (2006) Charles P. A., Coe M. J., 2006, in Compact Stellar X-ray Sources (eds Lewin W. H. G. & van der Klis, M.), Cambridge Univ. Press, Cambridge, UK, 215–265
  * Cherepashchuk et al. (2005) Cherepashchuk A.M. et al., 2005, A&A, 437, 561
  * Fabrika (2004) Fabrika S., 2004, Astrophys. Space Phys. Rev., 12, 1
  * Falcke & Biermann (1995) Falcke H., Biermann P., 1995, A&A, 293, 665
  * Gaisser (1990) Gaisser T. K., 1990, Cosmic Rays and Particle Physics, Cambridge University Press, Cambridge
  * Ghisellini, Maraschi & Treves (1985) Ghisellini G., Maraschi L., Treves A., 1985, A&A, 146, 204
  * Gies et al. (2008) Gies D. R. et al., 2008, ApJ, 678, 1237
  * Ginzburg & Syrovatskii (1964) Ginzburg V. L., Syrovatskii S. I., 1964, The origin of cosmic rays, Pergamon Press Ltd.
  * Kelner, Aharonian & Bugayov (2006) Kelner S. R., Aharonian F. A., Bugayov V. V., 2006 Phys. Rev. D, 74, 034018
  * Kelner, Aharonian & Bugayov (2009) Kelner S.R., Aharonian F.A., Bugayov V. V., 2009, Phys. Rev. D, 79, 039901
  * Khangulyan et al. (2007) Khangulyan D. et al., 2007, MNRAS, 380, 320
  * Körding, Fender & Migliari (2006) Körding E. G., Fender R. P., Migliari S., 2006, MNRAS, 369, 1451
  * Kosmas & Smponias (2018) Kosmas O. T., Smponias T., 2018, Adv. High Energy Phys. 960296, arXiv:1808.00303 [astro-ph.HE].
  * Kosmas (2020) Kosmas O. T., 2020, in preparation
  * Lipari, Lusignoli & Meloni (2009) Lipari P., Lusignoli M., Meloni D., 2009, Phys. Rev. D, 75, 123005
  * Long et al. (1981) Long K. S., Dodorico S., Charles P. A., Dopita M. A., 1981, ApJ, 246, L61
  * Mirabel & Rodríguez (1999) Mirabel I.F., Rodríguez L.F., 1999, ARA&A,37, 409
  * Orosz (2003) Orosz J. A., 2003, in A Massive Star Odyssey: From Main Sequence to Supernova (eds van der Hucht K. A., Herrero A., Esteban C.), Proc. IAU Symp. 212, ASP, San Francisco, 365
  * Orosz et al. (2008) Orosz J.A., McClintock J.E., Narayan R. et al., 2007, Nature, 449, 872
  * Orosz et al. (2011) Orosz J.A. et al., 2011, ApJ, 742, 84
  * Papadopoulos (2020) Papadopoulos D. A., 2020, Study of gamma-ray and neutrino emission from microquasar jets, MSc Thesis, University of Ioannina (unpublished)
  * Papadopoulos, Papavasileiou & Kosmas (2020) Papadopoulos D.A., Papavasileiou Th.V., Kosmas T.S., 2020, J.Phys.Conf.Ser., to appear, [astro-ph.HE] arXiv:2010.00396.
  * Pietsch et al. (2006) Pietsch W. et al., 2006, ApJ, 646, 420
  * Reid et al. (2011) Reid M.J. et al., 2011, ApJ, 742, 83
  * Remillard & McClintock (2006) Remillard R. A., McClintock, 2006, ARA&A, 44, 49
  * Reynoso, Romero & Christiansen (2008) Reynoso M.M., Romero G.E., Christiansen H.R., 2008, MNRAS, 387, 1745
  * Reynoso & Romero (2009) Reynoso M.M., Romero G.E., 2009, A&A, 493, 1
  * Romero et al. (2003) Romero G.E., Torres D.F., Kaufman Bernadó M.M., Mirabel I.F., 2003, A&A, 410, L1
  * Romero & Vila (2008) Romero G. E., Vila G. S., 2008, A&A, 485, 623
  * Romero et al. (2016) Romero G. E., Boettcher M., Markoff S., Tavecchio F., 2016, Space Sci. Rev., 207, 5
  * Romney et al. (1987) Romney J.D. et al., 1987, ApJ, 321, 822
  * Saito et al. (2009) Saito T. Y. et al. (for the MAGIC Collaboration), 2009, preprint (arXiv:0907.1017v2)
  * Smponias & Kosmas (2011) Smponias T., Kosmas T. S., 2011, MNRAS, 412, 1320
  * Smponias & Kosmas (2014) Smponias T., Kosmas T. S., 2014, MNRAS, 438, 1014
  * Smponias & Kosmas (2015) Smponias T., Kosmas O. T., 2015, Adv. High Energy Phys., 921757
  * Smponias & Kosmas (2017) Smponias T., Kosmas O. T., 2017, Adv. High Energy Phys., 4962741, arXiv:1706.03087 [astro-ph.HE]
  * Spencer (1984) Spencer R. E., 1984, MNRAS, 209, 869
  * Tsoulos, Kosmas & Stavrou (2017) Tsoulos I. G., Kosmas O. T., Stavrou V. N., 2019, Comp. Phys. Commun. 236, 237
  * Vieyro & Romero (2017) Vieyro F. L., Romero G.E., 2012, A&A, 542, A7, 1

##  Appendix A

###  A.1 The hot proton’s injection function

At the observer’s frame of reference, the injection function of protons,
$Q_{p}(E,z)$, and is given by (Reynoso, Romero & Christiansen, 2008):

$\displaystyle
Q_{p}(E,z)=\left(\frac{z_{0}}{z}\right)^{3}\frac{Q_{0}}{\Gamma_{b}\left(E-\beta_{b}\cos\theta\sqrt{E^{2}-m^{2}c^{4}}\right)^{2}}$
$\displaystyle\times\left[1-\frac{\beta_{b}Ecos\theta}{\sqrt{E^{2}-m^{2}c^{4}}}\right]$
(36)

where $\Gamma_{b}$ is the bulk Lorentz factor of the jet, and $\theta$ denotes
the angle between the jet’s injection axis and the direction of the line of
sight (LOS).

Then, the normalization constant $Q_{0}$ is obtained by determining the power
$L_{p}$ of the relativistic protons (Reynoso, Romero & Christiansen, 2008),
given by

$L_{p}=\int_{V}d^{3}r\int_{E_{p}^{min}}^{E_{p}^{max}}E_{p}Q_{p}(E_{p},z)dE_{p}$
(37)

After performing the latter integration, the normalization constant $Q_{0}$
(see text) reads (Reynoso, Romero & Christiansen, 2008)

$Q_{0}=\frac{2c}{z_{0}}K_{0},\quad\quad
K_{0}=\frac{4q_{rel}L_{k}}{cr_{0}^{2}ln\left({E_{p}^{\prime max}/E_{p}^{\prime
min}}\right)}\,.$ (38)

###  A.2 The distribution of pions produced per p-p collision

The distribution of pions produced per p-p collision is

$\displaystyle F_{\pi}^{(pp)}\left(x.\frac{E}{x}\right)=4\alpha
B_{\pi}x^{\alpha-1}\left(\frac{1-x^{\alpha}}{1+rx^{\alpha}(1-x^{\alpha})}\right)^{4}$
$\displaystyle\times\left(\frac{1}{1-x^{\alpha}}+\frac{r(1-2x^{\alpha})}{1+rx^{\alpha}(1-x^{\alpha})}\right)\left(1-\frac{m_{\pi}c^{2}}{xE_{p}}\right)^{\frac{1}{2}}$
(39)

with $x=E/E_{p}$, $B_{\pi}=\alpha^{\prime}+0.25$,
$\alpha^{\prime}=3.67+0.83L+0.075L^{2}$, $r=2.6/\sqrt{\alpha^{\prime}}$, and
$\alpha=0.98/\sqrt{\alpha^{\prime}}$ (Kelner, Aharonian & Bugayov, 2006,
2009).

###  A.3 Produced $\gamma$-ray spectrum

Following the treatment of (Kelner, Aharonian & Bugayov, 2006, 2009), the
spectrum of produced gamma-rays with energy $x={E_{\gamma}}/{E_{p}}$ for a
primary proton with energy $E_{p}$ is written as

$\displaystyle F_{\gamma}(x,E_{p})$ $\displaystyle=$ $\displaystyle
B_{\gamma}\frac{lnx}{x}\left[\frac{1-x^{\beta_{\gamma}}}{1+k_{\gamma}x^{\beta_{\gamma}}(1-x^{\beta_{\gamma}})}\right]^{4}$
$\displaystyle\times$
$\displaystyle\left[\frac{1}{lnx}-\frac{4\beta_{\gamma}x^{\beta_{\gamma}}}{1-x^{\beta_{\gamma}}}-\frac{4k_{\gamma}\beta_{\gamma}x^{\beta_{\gamma}}(1-2x^{\beta_{\gamma}})}{1+k_{\gamma}x^{\beta_{\gamma}}(1-x^{\beta_{\gamma}})}\right]$

where

$\displaystyle B_{\gamma}$ $\displaystyle=1.3+0.14L+0.011L^{2},$
$\displaystyle\beta_{\gamma}$ $\displaystyle=\frac{1}{0.008L^{2}+0.11L+1.79},$
$\displaystyle k_{\gamma}$ $\displaystyle=\frac{1}{0.014L^{2}+0.049L+0.801}$

with

$L=ln\left(E_{p}/(1TeV)\right)$ (42)