Abstract
The hadronic jets in a microquasar stellar system are modeled with the relativistic hydrocode PLUTO. We focus on neutrino emission from such jets produced by fast proton (nonthermal) collisions on thermal ones within the hadronic jet. We adopt a semianalytical 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 made the last decades. In contrast to the preset 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 hydrocode. 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.
1. Introduction
In binary stars commonly known as microquasars (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 multiwavelength and also particle emitter [2–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 microquasars (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–7].
In hadronic microquasar jets, the proton-proton interactions with the subsequent decay 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 () decay, as discussed in previous works [7–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 and mesons and neutrinos and 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 microquasars (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 (nonthermal) and cold (thermal) protons within the jets themselves [2–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, microquasar neutrino and gamma ray sources need to be modeled fully relativistically [1, 7]. A suitable treatment is offered by the relativistic hydrocodes developed recently, such as the relativistic magnetohydrodynamical (RMHD) PLUTO hydrocode [15] employed in [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 [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, for example, the IceCube [19] and KM3NeT [20], thus clarifying the potential for observing neutrino emissions from microquasars.
2. 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 occurring within the MQ jets themselves [2–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 ) 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 scale.
For the particle (proton) acceleration rate at shocks (first-order Fermi mechanism), we havewhere denotes the magnetic field and denotes the proton energy ( and are the usual parameters, i.e., the proton charge and the speed of light, resp.). The acceleration efficiency parameter in our present calculations is set equal to (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 () production channel, andfor the charged-pion () production channels, where , comprises and pairs, respectively.
Subsequently, the neutral pions and other neutral mesons decay quickly producing high-energy gamma rays. The charged pions (), needed for the purposes of our present work (and 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 Decay
From inelastic p-p scatterings among nonthermal protons and thermal ones within the hadronic jet, neutrinos are mainly produced through charged-pion decay (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 and a half-life of s) decay due to the weak interaction, the primary decay mode of which (with a probability of ) is a reaction leading to an antimuon (muon) and a muonic neutrino (muonic antineutrino) asA less important decay mode of (), with probability of occurrence just , is its decay into a positron (electron) and an electron neutrino (electron antineutrino) as In this work, we neglect the neutrino production through the latter channels.
2.1.2. Delayed Decay Channel (Delayed Neutrinos)
The other important source of neutrinos in hadronic jets is the decay mode of the produced muons (muon leptonic decay) in 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 [11] within a few percent over a large range of the fraction of the energy of the incident proton () transferred to the secondary particles, that is, the ratio , with 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 in studying stellar objects as neutrino and gamma ray sources; for example, the SS-433 system is widely known from the early 1980s 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 modeled 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, and so forth, 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.
3. Results and Discussion
The main results of this work refer to the mean number density of the nonthermal protons (obtained with the algorithms mentioned before and the PLUTO hydrocode), the pion injection function, and the pion energy distribution describing the pion governing (4). The evaluation of the emissivity of the prompt neutrinos relies on these calculations.
3.1. Nonthermal Proton Density
We begin our calculations by considering the production of nonthermal protons in the jet. The nonthermal 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 ultrarelativistic velocities. While in our previous studies we adopted a fast (nonthermal) proton jet density , equal to a tiny fraction () of the corresponding thermal proton density, in the present work, we assume a power-law distribution of the form , with [3]. In addition, we considered a spatial density distribution , coming out of explicit calculations with the PLUTO hydrocode 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 ), that is, the slow proton density , 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 snapshots which are plotted in Figure 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 (). 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 in 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 [4, 11], incorporating mainly synchrotron and adiabatic energy loss mechanisms.
3.2. Pion Injection Function and Pion Energy Distribution
For every p-p interaction (one “fast,” nonthermal 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; that is, we get a spectrum of possible energies for the resulting pion. That spectrum, per p-p collision, is represented by 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 [11], the function is given by the expression which represents the pion spectrum per proton-proton interaction. , , , , , and is the jet’s luminosity (see [4, 11]). In Figure 2(a), the product is plotted as a function of the ratio , for three different incoming fast proton energies ( GeV, GeV, and GeV), which cover the energy range of interest.
(a)
(b)
With the aid of this function, we calculate the pion injection function, , through the relationwhere . stands for the fast proton density, is the ratio of the pion energy to proton energy, and is the proton-proton inelastic collision cross section.
The pion injection function, , depends on the thermal proton density, . In Figure 2(b), we plot versus the pion energy for three different jet densities (, , and ). We notice the approximate square dependence of the scale of on the jet density, which is because also depends on .
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 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, that is, the pion injection function .
In order to obtain the pion distribution entering neutrino emissivity, we solve the following transport equation: where denotes the pion energy distribution. The numerical integration of the transport equation, for a cell of the hydrocode, that is, 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 , in the energy range of our interest, is illustrated in Figure 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 the stochastic portion of the general transport equation inactive. 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]).
3.3. Neutrino Emissivity
As mentioned before, in this work, we consider neutrinos emanating from direct pion decay (prompt neutrinos; see reaction (4)). In the semianalytical approach implemented in this work, the emissivity of prompt neutrinos is obtained with the aid of from the expression [4, 12] where and is the pion decay time scale. is the well-known theta function (for further parameter details, see [7]). The neutrino emission calculation could be performed mainly following the analysis of [3, 4, 11, 12].
For the readers’ convenience, we should mention the following. The nonthermal proton distribution suffers synchrotron and adiabatic losses, affecting the balance in the transport between protons and pions. The total neutrino emissivity can then be calculated by adding up contributions from every volume element (3D cell) and dividing the sum by the area of a sphere with radius equal to 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 versus the neutrino energy, in Figure 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 Figure 4 are unnormalized, 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 (nonthermal) proton power in the jet during the eruption modeled. The latter energy is, in turn, a fraction of the total jet kinetic power, or kinetic luminosity, . For erg/s, then erg/s, if . Assuming now a neutrino fraction of , we then obtain erg/s.
In general, to scale the neutrino intensity, a normalization factor is needed. In our case, that factor results from energetic arguments. This factor should multiply the neutrino emissivity of Figure 4 to obtain the neutrino intensity along the jet. As an example, assuming, as above, a fast proton energy fraction of , by equating the area under the curve of Figure 4 ( corresponds to the average bulk proton density at the jet base) to the fast proton fraction of the jet kinetic luminosity ergs/s, we obtain approximately equal to erg/GeV. The latter corresponds to a maximum neutrino intensity of about , which is compatible with the results of [4] (see, e.g., Figure 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].
We should point out that, in order to convert quantities from the jet reference frame to our rest frame, the calculational procedure can be, for example, 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.26 c, the average flow estimated for the jet. In addition, the line-of-sight direction is assumed to be constant all over the jet, at an angle of 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. Such detailed calculations will be included in a future publication.
4. 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 magnetohydrodynamic (RMHD) code (hydrocode). More specifically, the PLUTO code was executed incorporating a toroidal magnetic field component in the jet, resulting in 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 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 semianalytic 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 semianalytical 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 modeling of radio and gamma ray emission and in this work with efficient estimation of neutrino emission events originating from microquasar jets. For the observation of such neutrino fluxes, current terrestrial detectors (e.g., IceCube at South Pole) are in operation.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.