Abstract
In this work, we simulate -rays created in the hadronic jets of the compact object in binary stellar systems known as microquasars. We utilize as the main computational tool the D relativistic magnetohydrodynamical code PLUTO combined with in-house derived codes. Our simulated experiments refer to the SS433 X-ray binary, a stellar system in which hadronic jets have been observed. We examine two new model configurations that employ hadron-based emission mechanisms. The simulations aim to explore the dependence of the -ray emissions on the dynamical as well as the radiative properties of the jet (hydrodynamic parameters of the mass-flow density, gas-pressure, temperature of the ejected matter, high energy proton population inside the jet plasma, etc.). The results of the two new scenarios of initial conditions for the microquasar stellar system studied are compared to those of previously considered scenarios.
1. Introduction
The emissions of -rays, neutrinos, etc. within the jets of microquasars (MQs) have recently gained great interest among researchers seeking to understand the structure properties and evolution of X-ray binary systems [1–3].
Special interests appeared on the -ray emission mechanisms inside the hadronic jets, as the photon-hadron interactions [4, 5] and the hadron-hadron interactions [6, 7] as well as the -ray absorption that help to deepen our knowledge on microquasars evolution [8].
On the other hand, the strong magnetic field in the jets may significantly affect the total internal -ray and neutrino emissions by tuning several processes determining the high energy proton population (synchrotron radio emission, etc.). Therefore, magnetic field effects should be appropriately incorporated and treated in jet models [9, 10].
Recently, neutrinos from galactic microquasars, even though not being detected so far, have been modelled and several simulations have been performed towards this aim. Such modelling may support future attempts to detect them (see, e.g., [7, 11]).
Invariably, the jets of microquasars as well as in general the astrophysical jets may be described as fluid flow emanating from the vicinity of the compact object. Such a microquasar system is the SS433 X-ray binary consisted of a donor (companion) star and a compact stellar object which emits relativistic jets in various wavelength bands. Currently it is the only microquasar observed with a definite hadronic content in its jets, as verified from observations of spectral lines [1, 2, 9].
Radiative transfer calculations may be performed at every point in the jet (for a range of frequencies/energies, at every location) [12], providing the relevant emission and absorption coefficients. In such cases, finally a line of sight integration may derive synthetic images of jet -ray emission, at the energy-window of interest [12, 13].
The relativistic treatment of jets takes into account various energy loss mechanisms that occur through several hadronic processes [4–7]. 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.
Many authors consider that the proton-proton (p-p) collisions between fast high energy protons (nonthermal protons) and bulk-flow slow (thermal) protons constitute the dominant cooling process of the high energy proton population of the jet. This mechanism explains the main part of the -rays and neutrinos produced in the binary SS433 system. The acceleration of thermal protons (diffusive first-order shock acceleration) occurs above a minimum threshold proton energy [7].
Assuming a Maxwellian energy distribution for the “slow” protons, only a tiny portion of the total bulk proton jet flow, i.e., the fastest of them, may undergo diffusive shock acceleration and may jump to the fast proton population. Hence, the fast protons constitute a small fraction of the total jet proton density which subsequently produce -rays, neutrinos, etc. In this work, we assume that this is the dominant mechanism generating high energy -rays in the SS433 microquasar jets.
For the sake of completeness, we mention that another rather important mechanism has been suggested based on the hadronic interactions occurring within the jet-wind interaction zone [7]. In this scenario, the -rays are generated from the decay of neutral pions, as . Pions are created via inelastic collisions of jet protons, ejected from the compact object, and ions of the stellar wind (such a process may also occur in the vicinity of the extended disk of the binary system) [6]. The latter emission mechanism is rather weak in SS433 [7].
So far microquasar -ray emissions have been observed through Cherenkov telescopes (HESS, MAGIC, and CTA) and orbital telescopes (INTEGRAL, Fermi) [15–20]. We also mention that, for low energy -rays, ongoing and future or next generation measurements with INTEGRAL (ESA satellite) and Fermi (NASA orbital telescope) may provide new data. Furthermore, very high energy -rays, in general above about GeV, can be studied with ground-based Cherenkov telescopes [21].
Phenomenologically, estimations of high energy -ray emission from MQs have extensively been carried out [9, 12]. In this work, using the D relativistic hydrocode PLUTO [22] and some in-house (mainly radiative transfer) code (now written both in Mathematica and in C) [23–25], we model -ray emissions from hadronic microquasar jets in the -energy range TeV.
The emission/absorption coefficients are computed on the basis of Monte Carlo simulations of terrestrial particle-particle collision experimental data [12, 26, 27] that describe -ray emission in MQs. Such simulations provide analytical parametrization for emission and absorption coefficients in a wide range of -ray energies (frequencies) produced in microquasar jets [12, 28].
Furthermore, by exploiting the hydrodynamic variable values supplied by PLUTO, our line of sight code may provide emission/absorption coefficients for every location in the jet. The results produced this way depend on the initial high energy proton distribution inserted in the hydrodynamical model jet [12, 29].
In the rest of the paper, at first (Section 2) the main MQs emissions mechanisms are briefly summarized. Then (Section 2.4), the radiative transfer method and the calculational procedure for obtaining gamma-ray emission are briefly described. The results of the 3D relativistic hydrocode PLUTO for the emission/absorption coefficients are presented in Section 4. Finally (Section 5), the main conclusions extracted in this work are summarized.
2. Outline of MQ Jet Emission Mechanisms
The SS433 microquasar, an eclipsing X-ray binary system with a compact object most likely a black hole, comprises two oppositely directed precessing hadronic jets. The spectrum of the companion (donor) star suggests that it is rather a late A-type MQ. In modelling -ray emission from SS433 in our present work, we assume that they are created mainly through p-p interactions between fast (relativistic) and slow (cold) protons within its hadronic jets.
Other production mechanisms, though not excluded, are considered less important. For example, some authors considered that the high energy -rays in hadronic MQ jets, are produced from p-p collisions taking place in the jets due to the interaction of relativistic protons with target protons of the rather weak stellar wind created in the companion star [6].
The main reaction chain that produces -rays starting from p-p interaction through the pion decay channel is written as( g and g). References [6, 7, 26] present an analytical description of the evolution of reaction chain within the jet. Here we assume that a very energetic but small proton population, (formed due to shock fronts in the jet), interacts with the bulk-flow jet protons.
From the latter protons, high energy protons are produced through first-order Fermi acceleration that occurs at shocks within the jet [7]. Such shocks are considered rather homogeneously distributed throughout the jet. The jet matter density is closely related to the density of the aforementioned shocks; thus, the internal shocks convert a portion of the bulk kinetic energy of cold protons, K, to the fast protons energy of the multidirectional motion. The rate, , at which some slow protons are transferred to the high energy distribution is described in [1] bywhere denotes the proton charge, with being the jet matter’s local velocity, and denotes the magnetic field.
Concerning the magnetic field , we assume that this is either constant or it decreases with the distance from the jet base as [30, 31]. We stress that such a variation leads to a decrease of proton acceleration not more than two orders of magnitude compared to its value around the jet’s base . We further stress that the acceleration rate of fast protons, , depends on the magnetic field as indicated in (2).
Alternatively, (2) gives the production rate of fast protons at every “location” in the jet, though the production of -rays from these fast protons occurs at a next stage (by “location” we mean a hydrodynamical grid-cell which, microscopically, is very large); of the order of cm [7]. We note that the presence of in (2) cuts off -ray emission from slow-moving matter into the jet (acceleration sites are much less in slow matter).
Moreover, in jet emission calculations the is incorporated into the jet density. Also, the proton acceleration rate cannot affect the -ray emission rate, unless drops below some value. As we will see below, in this work, instead of the proton density , we also use the product of jet matter density times the velocity squared [7, 13].
Regarding the particles ejected from a hadronic jet, we consider that they are mostly slow (thermal) protons of density and a small portion of fast (nonthermal) protons of density . The energy distribution of the fast protons in the jet’s frame (see below) is described by [7, 13]() which is a power law type distribution. The parameter takes the value and denotes a normalization constant [13].
In Table 1, we tabulate the values of some model jet parameters (together with explanation of their symbols) relevant to -ray emissions from the SS433 binary system (The scenarios C and D are described below).
2.1. Flux in Observation and Jet’s Frame
In our -ray flux calculations, we denote the flux density of the fast (slow) proton populations as (), in the observation frame, and as (), in the jet’s frame. Moreover, at a given jet point, in the jet (moving) frame, the fast protons energy-spectrum is described by (3). For the corresponding fast protons spatial density we adopt the relationwhere(, with being the magnitude of the total velocity vector) where is in protons/cm3 [29]. From (5) one can conclude that the emission from the fast-moving matter of the jet is larger compared to the emission from slow-moving matter which is because in (2), is proportional to and creates fast proton jet density which subsequently allows for p-p collisions to occur and for -rays to be produced. [6].
2.2. The Model of Jet’s Dynamics
The jet is assumed to travel along the axis (we consider particles of mass ). Then, for the flux densities, we can write (steady state)In general, is not necessarily parallel to the -axis, but it may point almost anywhere which in turn means that emission may occur from jet matter moving in any direction. Furthermore, the emission mechanism is based on “randomly oriented turbulent shocks”, so the emission is considered multidirectional (no secondary emissions from scattering are assumed, since more shocks exist wherever the jet matter moves faster).
In our simulations, large turbulences of the jet flow may appear which favor shocks existence. This is due to the assumed strong dependence on the local velocity of the jet or ambient matter (acceleration rate is proportional to and further ). Here, instead of the simple dependence, we adopt, in addition, the dependence to distinguish the moving matter of the jet from that of the surrounding medium. This way, the calculation of -rays and neutrino emissions from the jet are decoupled from the influence of the surrounding matter. Then, the jet’s contribution to high energy -ray emission is mostly dependent on its internal turbulence (turbulence here means spatial number density of proton accelerating shocks randomly oriented) [7].
From the above discussion, we note in short that the proton acceleration efficiency of the model jet is proportional to the square of the local velocity of the flow (the fast proton density is considered proportional to the square of the local velocity). Thus, the fast protons spatial density, , is also taken as proportional to the slow protons spatial density, as well as to the square of the local velocity.
Furthermore, for hydrodynamical jets [12, 13] the fast proton current density, , as a function of their energy, is given byIn the latter expression, denotes the slow bulk jet protons local density hydrodynamical model (PLUTO code) [12].
2.3. The Current Density in the Observer’s Frame
Regarding the transformation, to the observer frame, we write [32] represents a function of stationary frame energy written asIn the latter equation, denotes the angle between the jet axis and the line of sight (for SS433 microquasar ), andis the jet Lorentz factor.
Thus, provides the relation of (for laboratory frame) that depends on the -ray energy as measured in laboratory frame (see (7)). In conclusion, one can work with laboratory frame quantities only, which are also the jet model quantities (for other symbols the reader is referred to [12, 13, 28, 29]
2.4. The 3D Radiative Transfer in Time-Dependent Jet
The propagation of -rays along a one-dimensional line of sight (without scattering) we address here is based on the relationwhere is the absorption coefficient at a given frequency (energy) and is the relevant emission coefficient. denotes the intensity and is the length along the line of sight.
By considering the model jet artificially imaged in -rays, we calculate the emission from a small jet element corresponding to a computational cell (for simulated emissions). To this aim, we first define the quantity to represent the intensity created from a cell of volume , at a given frequency, or -ray energy , while is the hydrodensity of the cell () and stands for the emission coefficient of the cell at the same frequency. An alternative version of the above quantity iswhere is the local jet matter velocity.
3. Use of PLUTO Code for Gamma-Ray Emission Calculations
Our calculation of the emission coefficient proceeds directly starting from the hydrodynamical properties of the model jet (supplied by a check point of the PLUTO code). These quantities enter the calculation of the emission coefficients, , at every computational cell of the D hydrodynamical model grid.
We mention that the production of synthetic images from the data is carried out by using the line of sight code constructed in [12] (here, instead of the radio emission and absorption coefficients we require their -ray equivalents, and ; in the case of neutrino production we need only emission coefficients).
In using the hydrodynamical code, PLUTO, the energy (in GeV) refers to the observed -ray energy. The quantity , i.e., the bulk-flow slow jet proton number density of ejected particles, is the dynamically important. This number density is taken to represent the hydrodynamic number density of the PLUTO code, namely,The fast proton density, , though not-important dynamically, is radiatively important. For SS433, in the fast proton power law energy distribution, the index takes the value , and the ratio of the initial jet beam speed divided by the speed of light is [13].
In the hydrodynamic (HD) simulations with PLUTO code, as we have done previously [12, 13], the magnetic field lines are assumed to follow the matter flow. Their tangling with the jet material makes applicable the fluid approximation within the jet [33]. In this case, the magnetic field could not affect the flow dynamics which is however possible in the magnetohydrodynamic (MHD) treatment of PLUTO. In the relativistic hydrodynamical version of PLUTO, the magnetic field within the jet’s medium is assumed rather strong so as the coupling effects permit the fluid approximation to be applicable. At the same time, the dynamical effects of the magnetic field on the relativistic flow are not permitted [33].
We mention that, in modelling the microquasar SS433 system with PLUTO, only one of the twin jets is considered. The counterjet is presumed to exist outside the model space (at the bottom of plane), but its interference with our model system is considered very small [34]. The computational grid is D Cartesian , homogeneous and the boundary conditions are adopted to be reflective at the jet’s base ( plane) and outflow at all other planes of the computational domain (box).
The grid spans (for x, y, z, respectively), in model length units (equal to cm) and the resolution used is (for x, y, z, respectively. The jet emanates from the middle of the plane, at the point (60, 0, 60)cm and then advances while precessing around the (60, y, 60) line (parallel to the y-axis). The precession angle for the SS433 jet used ( radians), is slightly smaller than the value of degrees of [34], in order to allow the use of finer resolution.
The centre of the companion star is supposed to be outside of the box, at the point (400, 0, 400) while compact object is situated at the point (60, 0, 60), i.e., at the jet’s base. We remind that, because the exact orbital separation in SS433 is not well known, this estimation is within an order of magnitude (the objects are orbiting around their centre-of-mass). We also mention that we assume that the companion star is not included in the model [12]. However, its wind is included through its density which is taken as decreasing with distance r (as ), away from its centre.
Furthermore, we also include a simplified accretion disk wind through the jet’s dynamic interaction with both winds. This means that our model is less realistic as we approach the companion star and accretion disk locations but the results are reliable in the vicinity of the jet. For a detailed discussion related to important phenomena of the jet’s interaction zone with nearby winds, the reader is referred to the Refs. [12, 13, 28, 29] and references therein.
4. Results of Simulations for the New Jet Model Scenarios
In this section, we present the results for two new scenarios of initial conditions (referred to the microquasar stellar system SS433) obtained as follows: (i) Hydrodynamical simulation carried out by utilizing as main computational tool the 3D relativistic hydrocode PLUTO and (ii) Gamma-ray emission synthetic images obtained with the line of sight integration.
In scenario C, the jet precesses faster than reality, therefore precession effects are enhanced. In this case, the jet involves artificially accelerated precession, in order to better investigate the effects of precession on the surrounding winds, within the limited time-frame of the model run.
In scenario D, the jet is quite heavier than both winds, in order to consider the possibility of a dense jet beam, containing the estimated jet mass flow of SS433 while remaining more focused and more both narrow. The heavier jet of this scenario crosses the winds with greater ease. Also the effects of its interaction with the winds appear decreased.
We note that another characteristic scenario would have been a jet much lighter than both winds, but this would have taken longer simulation time, and practically more difficult.
In both cases, the jet begins to expand into the accretion disk wind, but at a more limited pace, due to the increased density of that wind in the model. As soon as the jet head reaches the stellar wind region, however, the jet’s expansion rate increases greatly (especially sideways), in the form of a side shock that accumulates ambient matter. At the same time, the accretion disk matter is expelled outwards, from the vicinity of the jet base, forming a “ring” around the jet. The accretion disk wind is swept in a prominent way, being denser than the stellar wind, leading to the creation of a halo around the jet base (see below).
The structure develops throughout the model run, therefore suggesting the possibility of its persistence later on, when the jet reaches its lobe in the W50 nebula. This is similar (to a certain extent) in structure to that discussed in [14]. The above scenarios are applied to the SS433 microquasar as described below.
4.1. Description of Runs for Scenarios C and D
4.1.1. Simulations of Scenario C
Scenario C (medium resolution, see Figures 1, 2, 3, and 4) has been chosen to cover the case of artificially fast precession of the jet. This way we may investigate the effects of precession on the system observables. The precessing jet sweeps across more of the ambient matter in a given period of time, as compared to an otherwise same but nonprecessing jet (run2-scenario), with both jet examples considered to be moving through identical surroundings. Consequently, the effects of the precessing jet on its surrounding environment are, in terms of affected volume, more prominent than when precession is absent. More ambient matter is displaced and part of it ends up being dragged along by the jet, albeit at a pace clearly slower than when precession is slower.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
The precessing jet model was therefore run, at an accelerated precession rate, and that showed an enhanced effect of sweeping the accretion disk wind matter from the jet cone. This, in effect, caused the formation of an enhanced, outward moving, “halo”, around the jet base (Figures 1 and 2). The precessing jet advances, first through the accretion disk wind, and then through the stellar wind, opening its path at an accelerating pace, due to meeting with progressively lower resistance, due to the falling density of the winds. The latter originated from both the accretion disk and the companion star. The gradient of the stellar wind, within the computational grid, is, however, not big, due to the increased distance from its origin, the companion star.
The jet precession pronounces the sweeping of the inner, denser part of the accretion disk wind, leading, later on, to the formation of an expanding approximately torus-shaped halo surrounding the jet whose axis roughly coincides with the axis around which the jet precesses. The torus consists mainly of accretion disk wind matter, forming a loose “barrel”, or hollow cylinder (Figures 1 and 2) around the jet, having a velocity component (clearly slower than the jet, i.e., subrelativistic) parallel to the jet’s precession axis, and a sideways expansion velocity component as well.
In the scale of the simulation, the companion star has a nonnegligible distance from the jet base, therefore in the immediate vicinity of the jet base, it is the more localized accretion disk wind that dominates over the stellar wind, in terms of density. It is, therefore, the accretion disk wind matter that is mainly expelled from the incoming jet and rushed outwards, swept over by the precessing jet. This finding might suggest a behaviour that over a much longer timescale leads to the formation of a “ruff” of material, around the cone swept by the precessing jet, perhaps along the lines of the “bow-tie” structure recently observed in SS433 [14].
The jet precession makes it possible to drag an increased quantity of surrounding matter along the jet which is in contrast to a nonprecessing jet where the swept matter is significantly less. This happens because the precessing jet covers a cone with an opening angle much larger than the jet’s own. Therefore, more ambient matter is displaced than from a straight jet. Furthermore, the partially sideways motion of the precessing jet further disturbs the surrounding winds, pushing and dragging them in an outwards direction. However, the time scale for precession is much bigger than the jet crossing time of the model space. Therefore, a longer term simulation, perhaps including replenishment of the winds as well, would be needed in order to study the effects of precession on the jet’s environment.
4.1.2. Simulations of Scenario D
In scenario D (medium resolution) discussed in this work (see Figures 5, 6, 7, and 8), the jet is assumed heavier than its surrounding winds, which in turn are also somewhat heavier than those of the other cases, leading to a faster crossing of the computational domain (Figure 5). Sideways expansion is also swift as the increased jet mass density allows for a faster “sweep” of the wind matter. Leftovers from the displaced accretion disk wind matter can be seen piling up around the jet base (Figure 6).
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
(d)
The jet forms a funnel that transfers mass outwards at an increased flow rate. The properties of the jet’s surroundings are now less pronounced as the expansion meets with reduced resistance from ambient matter. The dynamic behaviour of the jet dominates the hydrosimulation, with the wind’s matter giving way to the jet (Figures 7 and 8). This case covers the possibility of jet production as a very dense inflow at the source, thus enriching nearby interstellar matter with important mass outflow per time interval.
All of the hydrocode snapshots, of hydrodynamic parameters, have been created with the VisIt visualization code (https://wci.llnl.gov/simulation/computer-codes/visit), whereas the synthetic -ray images have been produced using IDL (https://www.harrisgeospatial.com/SoftwareTechnology/IDL.aspx).
4.2. Gamma-Ray Emission Synthetic Images
The above discussed runs with the hydrocode PLUTO have been performed for a precessing jet model of SS433. As the simulation proceeds, at some point the computational model space data is transferred to an output file to be processed (with the available line of sight (LOS) code) for producing a synthetic -ray image of the system. The data of a snapshot from the PLUTO hydrocode is transferred, in the form of 3D data arrays (density, velocity vector, pressure: ) to a routine that performs the LOS integration [12, 13]. Along the LOS, -ray emission and absorption coefficients are provided at each point.
At this level, the orbital separation of the compact object and the donor star is taken to be about cm so the stellar wind origin does not coincide with the jet base. Also, the jet is taken to travel through a halo produced from the accretion disk (centred at the jet’s base). This matter is commonly called the accretion disk wind [12, 13].
With our method we follow two separate steps of calculations. We first obtain the hydrodynamic quantities of density and velocity with the PLUTO hydrocode and second, the integration (with the LOS code) provides the intensity. This decoupling allows, in the calculation of emission, the LOS integration to be performed using either (in CGS units) or (in speed of light dimensionless units, ) as the emission coefficient. Then, the synthetic image is formed, and subsequently the values of all of its pixels are summed to provide the total intensity released from the studied object. Finally, the -ray emission calculation is performed separately in Mathematica (for unit proton number density).
In each of the above scenarios (runs), for -rays a synthetic image was produced, at a suitable model time, using the relevant radiative transfer code of [12] (only the emission coefficient was employed). As can be seen, in all the synthetic -ray images (Figures 4 and 8), at each point of the computational volume, the denser the matter, the higher the emission at that point is.
Furthermore, it is clear that the larger the number of significant emission points along a line of sight, the higher the total emission of the whole line of sight is. By adding the dependence on the local velocity, denser but slower matter cannot emit any significant amount of -rays. Therefore, emission from the jet body and its interaction zone with surrounding media can be seen to be stronger in the maps. In addition, the rest of the (roughly inert) medium in the system contributes very little to -ray emission.
Before closing we note that currently we apply the above-mentioned method to carry out jet emission simulations for other microquasar systems like the Cygnus X-1 and Cygnus . From the viewpoint of observations we should mention that, for lower energy -rays, orbital platforms, such as NASA’s Fermi and ESA’s INTEGRAL, already offered an important relevant body of observations for various systems [21, 35]. Very high energy -rays can also be studied using data provided by ground-based Cherenkov telescopes such as HESS, MAGIC and HEGRA. Hydrodynamical jet models, combined with artificial imaging, offer realistic estimates of the conditions in the jet and surrounding environments, adding insight to open questions about the -ray jet emission from a microquasar.
5. Summary and Conclusions
A precessing jet was modelled using relativistic hydrodynamic code (PLUTO). Furthermore, the results were processed using line of sight code, assuming that the flow velocity is much smaller than the speed of light. The LOS code integrates along lines of sight the equation of radiative transfer without scattering. The emission coefficients are, in general, a function of hydrodynamical and radiative parameters.
The intensity result of each LOS is assigned to the pixel where the LOS meets the imaging plane of “observation”. In this way, an image is formed which could be called a synthetic -ray image. Only emission is used for this paper, but absorption may both be incorporated.
The currently available resolution for -ray imaging is lower than the resolution of synthetic imaging. Yet the connections between the emission properties and the underlying system dynamics do offer useful constrains on a variety of system parameters, such as the jet kinetic luminosity energy . The above occurs in the light of potential future observations, originating from orbital -ray telescopes, from terrestrial Cherenkov detector arrays and even from underground neutrino detectors.
Data Availability
All data are provided in full in the results section of this paper.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
Dr. Odysseas Kosmas wishes to acknowledge the support of EPSRC via Grant EP/N026136/1 “Geometric Mechanics of Solids”.