Neutrino emission from magnetized micro-quasar jets

The hadronic jets in a micro-quasar stellar system are modelled with the relativistic hydrocode PLUTO. We focus on neutrino emission from such jets produced by fast proton (non-thermal) collisions on thermal ones within the hadronic jet. We adopt a semi-analytical approximation for the description of the secondary particles produced from p-p collisions and develop appropriate algorithms using the aforementioned injected protons as input. As a concrete example we consider the SS-433 X-ray binary system for which several observations have been performed the last decades. In contrast to the pre-set distribution of the fast protons along the jet employed in our previous works, in the present paper we simulated it by using a power-law fast proton distribution along the PLUTO hydro-code. This distribution gradually sweeps aside the surrounding winds, during the jet advance through the computational grid. As a first step, in the present work the neutrino energy spectrum is extracted from the model jet, facilitating a range of potential dynamical simulations in currently interesting microquasar jet systems.


Introduction
In binary stars commonly known as micro-quasars (MQs), two oppositely emitted jets of matter and radiation are produced. These systems are similar to Active Galactic Nuclei (AGN or quasars) and consist of a main sequence star (the giant companion or donor star), in coupled orbit with a compact astrophysical object (a neutron star or a black hole) [1]. A characteristic mass accretion disk develops close to the compact object from mass absorption through the inner Lagrangian Point (Roche Lobe Overflow) due to angular momentum conservation. The jets of a MQ appear quite collimated (due to the presence of a rather strong magnetic field) forming a multi-wavelength and also particle emitter [2,3,4]. Stellar MQs are currently important astrophysical systems with growing interest in their investigations within astrophysics, particle physics and cosmology. In the case of black hole micro-quasars (when the compact object is a black hole) the stellar system provides excellent testing grounds for black hole theories. Therefore, an improved understanding of the dynamical astrophysical conditions within the jets in MQs is of significant importance [5,6,7].
In hadronic micro-quasar jets, the proton-proton interactions with the subsequent decays of the secondary particles, mostly π ± mesons, produce high energy neutrinos. These collisions result also in the production of high-energy gamma rays, through the neutral pion (π 0 ) decay, as discussed in previous works [7,8,9,10,11]. Recent simulations of high energy p-p interactions in terrestrial laboratories provide quite accurate energy distributions of secondary products in the high energy range (above 100 GeV) and determine parametric expressions of energy spectra for secondary particles like π 0 and π ± mesons, neutrinos but also for gamma rays and electrons produced in inelastic p-p collisions [11,12]. Such distributions may also be implemented when studying the hadronic MQs as neutrino and gamma ray sources [7].
Among the hadronic models proposed for the energy emission from micro-quasars (MQs), two are the most important: (i) In the first, relativistic protons in the jet interact with target protons from the stellar wind of the companion star. (ii) In the second, neutrinos and gamma-rays are produced from p-p interactions between relativistic (non-thermal) and cold (thermal) protons within the jets themselves [2,3,4,11,12]. In the latter case, relativistic (fast) protons within the jet are subject to different mechanisms that can make them lose energy. It is interesting to know the energy range where p-p collisions are the main (dominant) cooling process that produces the corresponding neutrinos (or gamma-rays). On the other hand, the cold (slow) protons serve as targets for the relativistic protons [13,14].
From a phenomenological point of view, micro-quasar neutrino and gamma ray sources need to be modelled fully relativistically [1,7]. A suitable treatment is offered by the relativistic hydrocodes developed recently, such as the relativistic magnetohydrodynamical (RMHD) PLUTO hydro-code [15] employed in Ref. [7,16,17] in order to simulate the hadronic jets of the SS-433 MQ, an X-ray binary star [5,6,18].
The present paper, is an extension of our work of Ref. [7] where we modeled simulated neutrino emission from Galactic astrophysical hadronic jets originating from the vicinity of compact objects in binary stellar systems. Our dynamical simulations come out of the RMHD PLUTO code in conjunction with the in house developed (in C, Mathematica and IDL) codes. We now produce further results that aim to be directly comparable to the sensitivities of modern high energy neutrino detectors, e.g. the IceCube [19] and KM3NeT [20], thus clarifying the potential for observing neutrino emissions from micro-quasars.

Brief description of the main background and formalism
In this work, we adopt the model explaining the neutrino and gamma-ray production through the p-p interactions between relativistic and cold protons occuring within the MQ jets themselves [2,3,4,11,12]. Relativistic protons in the jet are subject to various mechanisms that can lead to energy release. As is well known, in the case of hadronic MQ jets, a small portion (about 1%) of the protons (bulk flow protons) may be accelerated through first order Fermi acceleration procedures that take place essentially at shock fronts inside the jet. In general, accelerated particles within the jet may gain energy up to the T eV scale.
For the particle (proton) acceleration rate at shocks (first order Fermi mechanism), we have where B denotes the magnetic field and E p the proton energy (e and c are the usual parameters, i.e. the proton charge and the speed of light, respectively). The acceleration efficiency parameter η in our present calculations is set equal to η = 0.1 (efficient accelerator case, mildly relativistic shocks near the jet base) [21]. From the scattering of fast protons off slow protons, high energy pions and kaons are produced which may further decay to very high energy gamma rays and neutrinos. The reaction schemes are described by equations of the form for the neutral-pion (π 0 ) production channel, and for the charged-pion (π ± ) production channels, where F i , i = 0, 1, 2, comprises π 0 and π + π − pairs, respectively. Subsequently, the neutral pions π 0 and other neutral mesons decay quickly producing high energy gamma-rays. The charged pions π + (π − ), needed for the purposes of our present work (also the charged-kaons), decay and lead to muons and furthermore to the production of various flavors of neutrinos as discussed below.
2.1. Secondary charged particle decays From inelastic p-p scatterings among non-thermal protons and thermal ones within the hadronic jet, neutrinos are mainly produced through charged pion decays (known as prompt neutrinos). The muons included in the by-products can afterwards decay again into an electron (or a positron) and the associated two light neutrino flavors (delayed neutrino beam) according to the reactions described below.
2.1.1. Prompt decay channels (prompt neutrinos): The π + (π − ) mesons (with a mass of m π = 139.6M eV /c 2 and a half-life of 2.610 −8 s) decay due to the weak interaction, the primary decay mode of which (with a probability of 0.999877) is a reaction leading to an anti-muon (muon) and a muonic neutrino (muonic anti-neutrino) as A less important decay mode of a π + (π − ), with probability of occurrence just 0.000123, is its decay into a positron (electron) and an electron neutrino (electron anti-neutrino) as In this work we neglect the neutrino production through the latter channels.

Delayed decay channel (delayed neutrinos):
The other important source of neutrinos in hadronic jets is the decay mode of the produced muons (muon leptonic decays) in the reactions (4), which produces also two neutrinos described by the processes In general, the analytical formulae suggested from laboratory p-p collisions resemble the simulated distributions extracted in Ref. [11] within a few percent over a large range of the fraction of the energy of the incident proton (E p ) transferred to the secondary particles, i.e. the ratio x = E i /E p , with E i being the energy of the secondary particle (e.g. pion).
From an experimental point of view, for astrophysical gamma rays and neutrinos extremely sensitive detection systems have been developed [8,9,20]. These detectors sparked a renewed interest for studying stellar objects as neutrino and gamma ray sources, e.g. the SS-433 system is widely known from the early 80's as the only MQ with a verified hadronic jet content. We mention, for example, that observations of iron lines in the spectrum of the SS-433 MQ provided useful information regarding the hadronic content of its jets [21].
From a theory and phenomenology point of view, the gamma ray and neutrino production from a hadronic MQ that are of interest in the present work, is based on reliably determining the distribution of the fast protons and the realistic injection functions of the produced secondary particles (pions, kaons, muons, etc.).
In previous works [7,16,17], the hadronic jet was modelled using the PLUTO code. The results of PLUTO were then processed in order to calculate the emissivity of various secondary particles (pions, kaons) and the produced muons, gamma-rays, etc., on the basis of the spatial and time variation of physical parameters like the magnetic field that collimated the jet, the mass number density for every grid cell of the PLUTO code and others.
Before proceeding to the presentation and discussion of the results we should mention that, the discrimination of prompt and delayed neutrinos from MQ jets is not possible, therefore the results obtained in the present work refer to physical quantities pertaining to prompt neutrinos, nonetheless as they are much faster to simulate computationally.

Results and discussion
The main results of this work refer to the mean number density of the non-thermal protons (obtained with the algorithms mentioned before and the PLUTO hydrocode), the pion injection function and the pion energy distribution describing the pion governing Eq. (4). The evaluation of the emissivity of the prompt neutrinos relies on these calculations.

Non-thermal proton density
We begin our calculations by considering the production of non-thermal protons in the jet. The non-thermal proton population emerges from the bulk jet flow that comprises mainly thermal protons, moving mildly relativistically. Some of the slow protons are locally accelerated, at shock fronts appearing within the jet flow (first-order Fermi acceleration process), to ultra-relativistic velocities. While in our previous studies we adopted a fast (non-thermal) proton jet density N p , equal to a tiny fraction (10 −6 ) of the corresponding thermal proton density, in the present work, we assume a power-law distribution of the form N p = N 0 E −α , with α ≈ 2 [3]. In addition, we considered a spatial density distribution n z , coming out of explicit calculations with the PLUTO hydro-code as discussed below.  For an RMHD simulation of a rather laterally restricted magnetized jet we, first, calculated the mean matter density along the jet axis (as a function of z), i.e. the slow proton density n(z), by evaluating the PLUTO density over a slice cut perpendicular to the jet axis. In order to cover the temporal evolution of the jet as the simulation evolves, these mean density values have been obtained for a number of 8 snapshots which are plotted in Fig. 1. From this Figure we can see how the mean density profile evolves along the jet. Its peak is moving outwards while the overall maximum gradually decreases. The jet remains confined, mainly due to the presence of a toroidal magnetic field component (B tor ). The surrounding wind helps shape the jet as well, especially at the early stages of the simulation, before the wind begins to be swept by the jet.
As the jet advances through the computational grid, it gradually sweeps aside the surrounding winds resulting to a near-steady-state with a rather flat density profile. The magnetic jet confinement prevents the jet density from falling too much along the jet. It is worth mentioning that, for the characteristic time scales of the energy loss mechanisms, we largely follow Refs. [11] and [4], incorporating mainly synchrotron and adiabatic energy loss mechanisms.

Pion injection function and pion energy distribution
For every p-p interaction (one 'fast', non-thermal proton scattered off a 'slow', thermal one) we obtain a probability density of a resulting pion at every position along the possible spectrum of resulting pions, i.e. we get a spectrum of possible energies for the resulting pion. That spectrum, per p-p collision, is represented by F π and is dependent on the incoming fast proton energy (slow proton Energy is negligible by comparison) and the ratio of a given position at the pion spectrum to the incoming proton energy. In Ref. [11] the function F π is given by the expression which represents the pion spectrum per proton-proton interaction. x = E/E p , B π = a + 0.25, a = 3.67 + 0.83L + 0.075L 2 , r = 2.6/ √ a , α = 0.98/ √ a , and L is the jet's luminocity (see [4,11]). In Fig. 2 (left panel) the product xF π is plotted as a function of the ratio x, for three different incoming fast proton energies (E p =10 3 GeV, E p =10 4 GeV, E p =10 5 GeV), which cover the energy range of interest.
With the aid of this function, we calculate the pion injection function, Q (pp) π through the relation where k = E/E (max) p . N p stands for the fast proton density, x is the ratio of the pion energy to proton energy, and σ inel pp is the proton-proton inelastic collision cross section. The pion injection function, Q (pp) π depends on the thermal proton density, n(z). In Fig. 2 (right panel), we plot the Q (pp) π versus the pion energy E π for three different jet densities (n=10 9 , n=10 10 , n=10 11 ). We notice the approximate square dependence of the scale of Q (pp) π on the jet density, which is because N p also depends on n(z).
As a physical interpretation, let us consider a large number of p-p collisions. So, we add up, at every pion spectrum energy, the contributions to the probability that a pion will result at that energy. Depending on the incoming proton energy for each collision, there may be a smaller or a larger contribution to any given pion energy, as long as it is smaller than the proton's , as a function of the ratio x=E π /E p . We consider charged pions (π ± ) as needed for our purposes in this work. E π denotes the secondary particle (pion) energy. Right: Variation of the pion injection function, Q π (E) through the pion energy spectrum, for three different jet densities n(z), located at different points along the model jet axis. energy in the first place (pion energy cannot exceed proton energy). So we integrate over many p-p collisions to find the pion spectrum of a collection of p-p collisions, i.e. the pion injection function Q (pp) π . In order to obtain the pion distribution entering neutrino emissivity, we solve the following transport equation where N π (E, z) denotes the pion energy distribution. The numerical integration of the transport equation, for a cell of the hydrocode, i.e. a localized position in space, is given by the following expression where We note here that, the physical conditions within a cell are taken to be constant and also that the macroscopic physical parameters (density, pressure, etc) within each cell are taken to be constant. Under these assumptions, the transport equation is only dependent on energy, which considerably simplifies its calculation. We also take the characteristic scale (mean free path) of the radiative interactions to be smaller than the cell size, leading to the containment of particle interactions within a given hydrocode cell. Furthermore, the time scale for the radiative interactions is so much smaller than the hydrocode's timestep, that the radiative interactions belong to a single timestep each time.
The behaviour of the pion distribution N π (E p ), in the energy range of our interest, is illustrated in Fig. 3. This curve refers to a typical computational cell of the PLUTO hydrocode. It could be easily extended to a number of hydrocode cells covering a span of the computational grid, therefore opening the way towards obtaining the neutrino emissivity from the whole grid.
In such a treatment, we consider a large number of interacting particles per computational cell, therefore the probability density in the transport equation can be approximated by the number density of the particles, rendering inactive the stochastic portion of the general transport equation. Moreover, only the deterministic portion of the transport equation is employed, which simplifies it to a deterministic partial differential equation (for further details on the meaning of various symbols and functions used in this section, the reader is referred to [4,11]).

Neutrino emissivity
As mentioned before, in this work we consider neutrinos emanating from direct pion decays (prompt neutrinos, see reaction 4). In the semi-analytical approach implemented in this work, the emissivity of prompt neutrinos, is obtained with the aid of N π (E p ) from the expression [4,12] where x = E/E π and t π,dec is the pion decay time-scale. Θ(χ) is the well-known theta function (for further parameter details see Ref. [7]). The neutrino emission calculation could be performed mainly following the analysis of Refs. [3,4,11,12]. For the readers' convenience, we should mention the following. The non-thermal proton distribution suffers synchrotron and adiabatic losses, affecting the balance in the transport between protons and pions. The neutrino emissivity can then be integrated over 3D cells volume (voxels volume) and divided by the surface of a sphere with radius the distance to Earth. The result is a synthetic 'neutrino emission observation' of the binary system. By repeating the process for many energies, we can then obtain a synthetic spectral emission distribution, for direct comparison with observations.
As an illustration of the behaviour of neutrino emissivity Q π→ν (E) versus the neutrino energy, in Fig. 4 the neutrino spectra from a series of computational slices, cut perpendicular to the jet axis (at equal intervals along the model jet), are shown. The density of each slice is spatially averaged over the slice surface and that average density (see Table 1) is then employed in the neutrino emission calculation. The averaging is performed in IDL and the emission calculation in Mathematica. The results presented in Fig. 4 are un-normalized, but in order to compare to minimum detection levels of existing and future instruments, the simulation results can be normalized energetically and then calibrated for a given specific instrument (this is going to be presented elsewhere).
The integrated (across the spectrum used) energy emitted, per unit time, through neutrinos, is presumed to be a fraction of the fast (non-thermal) proton power in the jet during the eruption modelled. The latter energy, is, in turn, a fraction of the total jet kinetic power, or kinetic luminocity, L k . For L k = 10 40 erg/s, then L f p = 10 36 erg/s, if f f p = 10 −4 . Assuming now a neutrino fraction of f n = 0.5 we then obtain L n = 0.5 × 10 36 erg/s.  In general, to scale the neutrino intensity, a normalization factor F norm is needed. In our case, that factor results from energetic arguments. This factor should multiply the neutrino emissivity of Fig. 4 to obtain the neutrino intensity along the jet. As an example, assuming, as above, a fast proton energy fraction of 10 −4 , by equating the area under the ρ 1 curve of Fig. 4 (ρ 1 corresponds to the average bulk proton density at the jet base) to the fast proton fraction of the jet kinetic luminosity L k = 10 40 ergs/s, we obtain F norm approximately equal to F norm =10 54 erg/GeV. The latter corresponds to a maximum neutrino intensity of about 0.5× 10 32 , which is compatible with the results of Ref. [4] (see e.g fig. 8 of this reference). We note that, even though our model is quite detailed dynamically, its level of detail cannot be fully compared to observations as of today. The reason is that current and upcoming terrestrial neutrino detectors cannot resolve that much detail, due to the distance of the microquasars from Earth. Therefore, when compared to observations, the predictive power of our model is not very different from that of simpler models [3,4].
Number density (in 10 10 proton/cm 3 ) We should point out that, in order to convert quantities from the jet reference frame to our rest frame, the calculational procedure can be e.g. that of [13] or that of [14]. In the present work, we apply the treatment of [13]. The jet direction has been incorporated as a global effect within the jet, by imposing a fixed angle between the velocity direction of the flow and the line of sight to an observer here on Earth. Furthermore, the jet flow speed is taken to be set to 0.26c, the average flow estimated for the jet. In addition, the line-of-sight direction is assumed constant all over the jet, at an angle of θ=78 degrees, to the jet axis. This was done to keep the calculations within limits. In principle, each computational cell may have a different setting for the angle between its local velocity and the line of sight, as well as for the local emission calculation performed using its localized velocity value. In both cases, a much longer computational time is required.
Before closing, it is worth mentioning that, the total power emitted from the jet obtained in the present work is only a first approximation to the intensity estimate. Consequently, a more detailed comparison to detectors is required. Our model is indeed able to provide, for a given direction to the observer, individual Doppler effects for each 3D computational cell, and then integrate them numerically (see Ref. [26]).

Summary and Conclusions
In the present work, we evaluated the emissivity of neutrinos originating from hadronic MQ jets, where p-p collisions occur at shock fronts, leading to cascades of secondary particles, culminating to neutrino emission. We have implemented a new model describing the mass distribution along the jet axis, using the PLUTO relativistic magneto-hydrodynamic (RMHD) code (hydrocode). More specifically, the PLUTO code was executed incorporating a toroidal magnetic field component in the jet, resulting to a confined jet structure, the degree of confinement depending on the value of the field. For each cross section slice, cut along the jet (perpendicular to the jet axis) we calculated the mean values of the mass density. Then, we proceeded, in this manner, to process a number of 100 slices, covering the spatial range from the jet base to the end of the computational grid.
The main conclusion extracted from this analysis was that the hydrocode model (not based on explicit geometrical assumptions), employed for the hadronic jet, is dynamically a realistic tool. This is why we decided to utilise the PLUTO code for dynamical calculations as a basis for further investigation of the neutrino and gamma-ray emissivities from the jet. For our present calculations, we used the semi-analytic approach, in order to estimate the neutrino emissivity, as described in our previous work.
In studying the neutrino emissivity per grid cell, we set up a model geometry reminiscent of the semi-analytical method, but using the PLUTO hydrocode, while employing the known radiative formalism as discussed in the Introduction. This computational tool has previously provided us with a realistic modelling of radio and gamma-ray emission and in this work with efficient estimation of neutrino emission events originating from micro-quasar jets. For the observation of such neutrino fluxes, current terrestrial detectors (e.g. IceCube at South Pole) are in operation.

Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.   Table 1) by integrating numerically Eq. (12). For the sake of comparison with observations and other predictions, the value of Q π→ν (E) should be multiplied by the normalization factor F norm =10 54 erg/GeV (see the text).