Simulations of neutrino and gamma-ray production from relativistic black-hole microquasar jets

Recently, microquasar jets have aroused the interest of many researchers focusing on the astrophysical plasma outflows and various jet ejections. In this work, we concentrate on the investigation of electromagnetic radiation and particle emissions from the jets of stellar black hole binary systems characterized by their hadronic content in their jets. Such emissions are reliably described within the context of the relativistic magneto-hydrodynamics. Our model calculations are based on the Fermi acceleration mechanism through which the primary particles (mainly protons) of the jet are accelerated. As a result, a small portion of thermal protons of the jet acquire relativistic energies, through shock-waves generated into the jet plasma. From the inelastic collisions of fast (non-thermal) protons with the thermal (cold) ones, secondary charged and neutral particles (pions, kaons, muons, $\eta$-particles, etc.) are created as well as electromagnetic radiation from the radio wavelength band, to X-rays and even to very high energy $\gamma$-ray emission. One of our main goals is, through the appropriate solution of the transport equation and taking into account the various mechanisms that cause energy losses to the particles, to study the secondary particle distributions within hadronic astrophysical jets. After testing our method on the Galactic MQs SS 433 and Cyg X-1, as a concrete extragalactic binary system, we examine the LMC X-1 located in the Large Magellanic Cloud, a satellite galaxy of our Milky Way Galaxy. It is worth mentioning that, for the companion O star (and its extended nebula structure) of the LMC X-1 system, new observations using spectroscopic data from VLT/UVES have been published few years ago.


I. INTRODUCTION
In recent years, astrophysical magnetohydrodynamical flows in Galactic and extragalactic microquasars (MQs) and X-ray binary systems (XRBs) have been modelled with the purpose of studying their multi-messenger emissions (e.g. neutrinos, gamma-ray emissions) [1][2][3]. For the detection of such emissions, extremely sensitive detector tools are in operation for recording their signals reaching the Earth like IceCube, ANTARES, KM3NeT [4][5][6], etc. Modelling offers good support for future attempts to detect radiative multiwavelength emission and particles, while in parallel several numerical simulations have been performed towards this aim [7][8][9].
Microquasars (or generally XRBs) are binary systems consisting of a compact stellar object (usually a black hole or a neutron star) and a donor (companion) star [10]. Mass from the companion star overflows through the system's Langrangian points to the equatorial region of the compact object gaining angular momentum and thus forming an accretion disc of very high temperature gas and matter. This mass expelling could be mainly due to a concentrated stellar wind [11][12][13]. Consecutively, a portion of the disc's matter is being concentrated by the system's magnetic field (initially it is attached to the ro- * Electronic address: th.papavasileiou@uowm.gr † Electronic address: odysseas.kosmas@manchester.ac.uk ‡ Electronic address: isinatkas@uowm.gr tating disc) and is ejected perpendicular to the disc in two opposite directions [14] forming the system's jets. The jets are detectable from the Earth even when the system's distance is too large. This is due to the relativistic velocities they acquire [15][16][17] combined with Doppler effects when they are headed towards the Earth. It has been shown that the kinetic luminosity coming from astrophysical microquasar jets constitutes a substantially large part of the total galactic cosmic radiation [18]. Moreover, it has been proved that microquasar binary systems that don't produce thermal jets are more likely to be neutrino and gamma-ray emission sources [19]. The most well-studied microquasar systems include the Galactic X-ray binaries SS433, Cyg X-1, Cyg X-3, etc. [20,21], while from the extragalactic systems we mention the LMC X-1, LMC X-3 (in the neighbouring galaxy of the Large Magellanic Cloud) [22,23], and the Messier X-7 (in the Messier 33 galaxy) [24]. Their respective relativistic jets are emission sources in various wavelength bands and high energy neutrinos. In this work, we focus our study on the extragalactic binary system LMC X-1 with the purpose to determine its gamma-ray and neutrino emissions produced through the processes and mechanisms that are about to be discussed. Concerning the binary system LMC X-1, new studies and observational results have been recorded, shedding more light on its unique spectral and environmental characteristics [25,26] and thus providing a more solid base for further in-depth analysis.
Up to now, the SS433 is the only microquasar observed with a definite hadronic content in its jets, as verified from observations of their spectra (see Ref. [7][8][9] and references therein). Radiative transfer calculations may be performed at every point in the jet for a range of frequencies (energies), at every location [27] providing the relevant emission and absorption coefficients. Line-of-sight integration, afterwards, provides synthetic images of γray emission, at the energy-window of interest [27,28]. Studies of the concentrations of jet particles, whose interactions lead to neutrino and gamma-ray production, need to take into account various energy loss mechanisms that occur due to several hadronic processes, particle decays and particle scattering [1,2]. In the known fluid approximation, macroscopically the jet matter behaves as a fluid collimated by the magnetic field. At a smaller scale, consideration of the kinematics of the jet plasma becomes necessary for treating shock acceleration effects.
In the model employed in this work [1,2], the jets are considered to be conic along the z-axis (ejection axis) with a radius r(z) = z tan ξ, where ξ its half-opening angle [29]. The jet radius at its base, r 0 , is given by r 0 = z 0 tan ξ, where z 0 is the distance of the jet's base to the central compact object. According to the jetaccretion speculation, only 10% of the system's Eddington luminosity [30] (L k = 1.2 × 10 37 M erg/s, M in solar masses) is transferred to the jet for acceleration and collimation through the magnetic field given by the equipartition of magnetic and kinetic energy density as B = 8πρ k (z) (see Ref. [1,2,27,28]).
By assuming the one-zone approximation [31], we consider a small portion of the hadrons (mainly protons) equal to q r ≈ 0.1 to be accelerated in a zone from z 0 to z max with the rate t −1 acc ≃ ηceB/E p [32] according to the 2nd order Fermi acceleration mechanism. The particles are accelerated to nearly relativistic velocities (with η = 0.1 being the acceleration efficiency) [33] resulting in a power-law distribution given in the jet's rest frame by In the rest of the paper, after a brief description (sec. II) of the p-p interaction chain leading to neutrino and gamma-ray emission, we solve the transfer equation (sec. III) to determine the energy distributions of the primary and secondary jet particles. Then (sec. IV), the obtained predictions for high-energy neutrino and gamma-ray production, for the LMC X-1 MQ system, are presented and discussed. Finally (sec. V), we summarize the main conclusions.

II. INTERACTION CHAIN LEADING TO NEUTRINO AND GAMMA-RAY PRODUCTION
In general, the main interactions of the relativistic protons include those with the stellar winds [34,35], the radiation fields (composed of internal and external emission sources) [36] as well as the cold hadronic matter of the jet. In this work, we focus on the last interaction because it is more dominant.
Neutrino production and gamma-ray emission are the results of a reaction chain that is caused by the p-p (and p-γ) interactions inside the jet. KM3NeT, ANTARES and IceCube [4][5][6] are prominent examples of undersea water and under-ice detectors that are able to detect those neutrinos that reach the Earth. The aforementioned reaction chain begins with inelastic p-p collisions of the relativistic protons with the cold ones inside the jet, which generate neutral (π 0 ) and charged (π ± ) pions according to the following reactions pp → pp + απ 0 + βπ + π − pp → pn + π + + απ 0 + βπ + π − pp → nn + 2π + + απ 0 + βπ + π − with α and β being the particle plurality factors that depend on the proton energy [37]. Neutral pions decay into gamma-ray photons as while the charged pions decay into muons and neutrinos as Subsequently, muons also decay into neutrinos. These are the main reactions feeding the neutrino and gamma-ray production channel in the employed model. All particles that take part in the neutrino and gammaray production processes lose energy while traveling along the acceleration zone, which could be due to different mechanisms as it is illustrated in Fig. 1. At first, the particles can be subjected to adiabatic energy losses due to jet expansion along the ejection axis with a rate depending on the jet's bulk velocity linearly [38]. Another important cooling mechanism, that depends on the cold proton density inside the jet, is due to inelastic collisions of the accelerated particles with the cold ones. The inelastic cross section for the p-p scattering is given by [7][8][9]39] where L = ln(E p /1000) with E p in GeV and E th = 1.2 GeV the threshold for the production of a single neutral pion. Respectively, for pions it holds [40] σ inel In addition, particles accelerated by the magnetic fields emit synchrotron radiation. Thus, they gradually lose Cooling rates for the relativistic protons and the secondary particles (π ± and µ ± ) produced after the p-p collision processes take place in LMC X-1 (first column) and Cygnus X-1 (second column).
part of their energy with a rate that heavily depends on the magnetic field and the particle energy.
Furthermore, protons and pions interact with radiation fields coming from internal or external regions of the jet, such as the system's accretion disc [41], as well as the synchrotron-emitting particles themselves. This leads to partially energy loss. Smaller particles, such as muons, transfer part of their energy to low-energy photons via inverse Compton scattering. However, such contributions can be ignored compared to those mentioned above.

III. SOLUTION OF THE TRANSFER EQUATION
The steady-state transfer equation, which fulfills the conditions discussed, is given by [7][8][9] where N (E, z) denotes the particle energy distribution in units of GeV −1 cm −3 while Q(E, z) is the particle source (or injection) function giving the respective production rate (in GeV −1 cm −3 s −1 ). Concerning the energy loss rate b(E), all the cooling mechanisms discussed are being represented as it is evident from its definition b(E) = dE/dt = −Et −1 loss . For each particle, the decay rate is added to the particle escape rate from the jet according to the relation: where z max stands for the end of the acceleration zone.
The general solution of the differential equation (4) is given by Therefore, in order to calculate the neutrino and gamma-ray emissivities, it is necessary to calculate first the energy distributions of all particles involved in the reaction chain (see sec. II).

Relativistic proton injection
In previous works [42], the appropriate injection function for the relativistic protons produced by the Fermi acceleration mechanism was found to be power-law with exponent ≈ 2 [43,44]. In the jet's rest frame, this powerlaw translates to the following expression where Q 0 is a normalization constant calculated through the total luminosity that is carried by the protons (or electrons) inside the jet (see Appendix A 1), where E min p = 1.2 GeV is the minimum proton energy that is sufficient and necessary for the Fermi mechanism to occur. The maximum energy is calculated by equating the particle acceleration rate with the total energy loss rate t −1 acc ≈ ηceB/E p = t −1 loss . For the binary systems of our interest, this is approximated as E max p ≃ 10 7 GeV . The dependence on z is due to particle conservation enforced on the respective current density [45]. Transformation of the source function of Eq. (6) to the observer's reference frame givens [8,9] Q(E, z) where Γ b responds to the jet's Lorentz factor and θ is the angle between the jet's ejection axis and the line of sight.

Pion energy distribution
The pion source function calculation requires the fast proton distribution as well as the p-p collision rate along with the pion spectra produced by each one of those collisions as where x = E/E p . In the latter integral, F π (x, E/x) denotes the pion mean number produced per p-p collision given by [39] F π x, where B π = α ′ + 0.25, α ′ = 3.67 + 0.83L + 0.075L 2 , r = 2.6/ √ α ′ and α = 0.98/ √ α ′ . Also, n(z) is the cold proton density of the jet written as where Γ is the cold proton Lorentz factor. From Eq. (8), it is worth noting that the proton distribution is entering both the neutrino and gamma-ray emissivity calculations through the pion injection rate Q π (E, z).

Muon spectra from pion decay
As it is known, for the muon energy distribution, both the mean right-handed and the mean left-handed muon numbers per pion decay are required for obtaining the total injection function. According to the CP invariance and provided that N π (E π , z) = N π + (E π , z)+N π − (E π , z), it holds [46]  where N + µ and N − µ represent the positive and negative right (or inversely the left) handed muon spectra, respectively (see Appendix A 2). In Eq. (11), we have used x = E µ /E π , r π = (m µ /m π ) 2 and Θ(y) the Heaviside function. Additionally, the pion decay rate is given by t −1 π,dec = (2.6 × 10 −8 γ π ) −1 s −1 , which implies that the pion distribution is important for the muon distribution calculation.

A. Neutrino spectra from pion and muon decay
After obtaining the energy spectra of the particles discussed above, one is able to estimate the total number of neutrinos produced directly from pion decay as well as from muon decay. Thus, the total emissivity contains both contributions as The first term represents the neutrino injection originating from pion decay as where x = E/E π , while the second term gives with x = E/E µ . In the latter equation, the muon decay rate depends on their energy as follows t −1 µ,dec = (2.2 × 10 −6 γ µ ) −1 s −1 . Also h 3 = h 4 = −h 1 = −h 2 = 1. From the four different integrals of the latter summation, the first and second represent the left-handed muons of positive and negative charge, respectively, while the third and fourth stand for the corresponding right-handed ones. Finally, integration over the acceleration zone gives the total neutrino intensity as [1,2] B. Gamma-ray emissivity for E > 100 GeV In this work, we assumed that gamma-ray production is mainly due to neutral pions decay which in turn are products of p-p inelastic collisions. The respective gamma-ray spectra have been simulated for photons of energy E γ = xE p to be [39] where B γ = 1.3 + 0.14L + 0.011L 2 , β γ = 1/(0.008L 2 + 0.11L + 1.79) and k γ = 1/(0.014L 2 + 0.049L + 0.801).
In addition, we have L = ln(E p /1T eV ). These results are consistent with proton energies in the range 0.1 T eV < E p < 10 5 T eV . Besides π 0 , Eq. (16) also considers the contribution of η mesons decay, which is approximately 25% when x ≈ 0.1.
In the energy range of our interest E γ ≥ 100 GeV , the gamma-ray emissivity, produced at a distance z from the compact object along the jet's ejection axis, is given by (in units of GeV −1 cm −3 s −1 ). The respective intensity is obtained again through integration over the acceleration zone. For E γ < 100 GeV, the delta-function approximation is employed.

C. Neutrino and gamma-ray intensity simulations
By introducing the calculated rates for the cooling mechanisms of particles that lead to the neutrino and gamma-ray production, to the transfer equation discussed previously along with the corresponding injection function, we are able to calculate the particle energy distributions. These results are shown in Fig. 2. The reaction chain seen above includes the relativistic protons accelerated by shock-waves through the Fermi mechanism, the charged pions produced by the p-p interactions and the muons that result from the pion decays. We consider all the above mechanisms and interactions taking place in the hadronic relativistic jets of the extragalactic binary system LMC X-1 located in our galaxy-neighbor LMC (Large Magellanic Cloud) and with a distance of 55 kpc to the Earth [22]. For comparison, the corresponding results for the galactic Cygnus X-1 are also shown. Then, the energy spectra of the produced high-energy neutrinos and gamma-rays can be simulated numerically by using those calculations.
The calculations were performed through the development of a C-code (that employs Gauss-Legendre numerical integration of the GSL library) and the use of the parameter values listed in Table I mostly describing geometric characteristics of the systems of interest.
At first, we calculated the fast proton distributions for three different distances to the central object for binary systems LMC X-1 and Cygnus X-1 as illustrated in the graphs a and b, respectively, of Fig. 2. It is evident that the total particle density production decreases while we move closer to the end of the acceleration zone even though the average particle energy increases.
In Fig. 2 graphs c and d, we show the energy distributions for pions that have suffered energy losses caused by  [49] various mechanisms such as synchrotron radiation emission, collisions with the rest of the jet matter etc. [42].
For comparison, we also plot the respective distributions of particles that do not lose energy at all. This exhibits the important effect which the cooling mechanisms cause on the total particle distributions. The same can be concluded for the muon distributions from the graphs e and f in Fig. 2. We also notice that the cooling rate characteristics are reflected upon the particle distributions as can be seen from the deviation of the two lines in the pion and muon case. Mainly, the deviation point coincides with the dominance of the synchrotron mechanism over the particle decay. The synchrotron losses take off causing the smoother transition in the solid line's case (with energy losses) compared to the dashed one (only decay) [42]. Furthermore, an instant observation would be that there is not significant difference between the two binary systems. This is a result of their black hole masses being of similar magnitude. The parameter that has the biggest impact though is the jet's half-opening angle, which roughly gives the degree of the jet's collimation. That is because the system's magnetic field, which is responsible for the collimation, presents a strong dependence on ξ.
After calculating all the necessary particle distributions, the neutrino and gamma-ray emissivities as well as the corresponding intensities are readily obtained in Fig. 3 and 4. For the LMC X-1 system, our results have shown that, the increase of the half-opening angle ξ leads to a decrease in the gamma-ray production, which is an expected result since the p-p collision rate drops with the jet's expansion as can be seen in Fig. 3.

V. SUMMARY AND CONCLUSIONS
During the last decades, the structure and evolution of relativistic astrophysical plasma outflows (jets) and specifically those connected to compact cosmic structures, became research subjects of intense interest. Towards this aim, a great research effort, experimental, theoretical, and phenomenological is absorbed by various high-energy phenomena, including production and terrestrial detection of high-energy cosmic radiation and neutrinos, originating from Galactic and extragalactic sources heading towards the Earth. We mention in particular, the jets of very large scale ejected from galaxy's quasar systems involving supermassive black holes at the central region (AGN) and those of much smaller scale involving stellar mass black holes and a companion star, known as micro-quasars and X-ray binary systems.
In this work, we concentrated on the latter class, the two-body systems consisting of a central object (usually a stellar mass black hole) and a companion star often of Otype or B-type main sequence stars. The former absorbs mass from the latter, forming an accretion disc of gas and matter which emits X-ray radiation due to very high temperatures prevailing in the area of accretion disk and black hole. The jets are formed when the system's magnetic field collects matter that is ejected away from the system in a collimated and accelerated bulk-like plasma flow.
The mass outflow can acquire relativistic velocities and is expelled perpendicular to the disc's surface. In many models these jets are treated magnetohydrodynamically by assuming various reliable approximations. In the present work, we considered the jet's matter to be mainly hadronic, with a portion of it accelerated through shock-waves to relativistic velocities. From the inelastic collision of relativistic protons on the cold ones (p-p interactions mechanism), secondary neutral and charged particles (pions, kaons, muons, etc.) are produced, the decay of which leads to high-energy neutrino and gammaray emissions.
One of our main goals was the calculation of the energy-spectra of high energy neutrino and gamma-ray produced inside such astrophysical jets. As concrete examples, we have chosen the Galactic X-ray binary Cygnus : Neutrino (a) and gamma-ray (b) intensities produced by pion (π ± ) decay that in turn are products of the p-p interaction mechanism taking place in the extragalactic binary system LMC X-1 as well as the well-studied system Cygnus X-1. X-1 system and the extragalactic LMC X-1 binary in order to simulate their neutrino and gamma-ray intensities emitted. For the observation of such high energy cosmic radiations and particle (neutrino) emissions, extremely sensitive detection instruments are operating and next generation detectors have been designed at the Earth like the IceCube detector (deep under the ice at South Pole), the ANTARES and KM3NeT (underwater in the Mediterranean sea), the CTA etc.
p /E min p ) (A1) The above result is calculated through the total luminocity carried by the protons that is given by

Right and left-handed muon spectra
The right-handed positive and negative muon spectra produced by pion decay are where x = E µ /E π and r π = (m µ /m π ) 2 . According to CP invariance, the number of µ − L produced by the π − decay is the same as the µ + R produced by the π + decay. Therefore, as the pion energy distribution refers to the sum of both π + and π − distributions, it is N + µ R = N − µ L . Similarly, it holds the same for the negative right-handed muons.