Next Article in Journal
Optimization of Quantum Noise in Space Gravitational-Wave Antenna DECIGO with Optical-Spring Quantum Locking Considering Mixture of Vacuum Fluctuations in Homodyne Detection
Previous Article in Journal
On the Possible Asymmetry in Gamma Rays from Andromeda Due to Inverse Compton Scattering of Star Light on Electrons from Dark Matter Annihilation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Simulated Radio and Neutrino Imaging of a Microquasar

by
Theodoros Smponias
Directorate of Primary Education of the Ionian Islands, 49100 Corfu City, Corfu, Greece
Galaxies 2023, 11(6), 110; https://doi.org/10.3390/galaxies11060110
Submission received: 31 July 2023 / Revised: 18 October 2023 / Accepted: 30 October 2023 / Published: 9 November 2023

Abstract

:
Microquasar stellar systems emit electromagnetic radiation and high-energy particles. Thanks to their location within our own galaxy, they can be observed in high detail. Still, many of their inner workings remain elusive; hence, simulations, as the link between observations and theory, are highly useful. In this paper, both high-energy particle and synchrotron radio emissions from simulated microquasar jets are calculated using special relativistic imaging. A finite ray speed imaging algorithm is employed on hydrodynamic simulation data, producing synthetic images seen from a stationary observer. A hydrodynamical model is integrated in the above emission models. Synthetic spectra and maps are then produced that can be compared to observations from detector arrays. As an application, the model synthetically observes microquasars during an episodic ejection at two different spatio-temporal scales: one on the neutrino emission region scale and the other on the synchrotron radio emission scale. The results are compared to the sensitivity of existing detectors.

1. Introduction

Microquasars (MQs) include a binary stellar system, with a main sequence star orbiting a collapsed stellar remnant [1]. Material from the star accretes onto a compact object, resulting in the production of a pair of relativistic ejections moving in opposite directions largely perpendicular to the binary orbital plane. These ejecta form jets, which emit anywhere between radio waves and very-high-energy (VHE) γ rays and neutrinos [2,3,4,5,6,7,8,9,10,11,12,13]. As shown in [2], the apparent superluminal motion in MQs suggests that jets comprise bulk hadron flows. The assumption of equipartition [6] leads to high magnetic field estimates, especially for inner jets [14]. The latter, together with the fluid approximation for the jet matter due to the presence of tangled magnetic fields [15,16], lays the foundation for the jet magnetohydrodynamic (MHD) approximation. A toroidal magnetic field may contribute to jet collimation along its path [14,17], while confinement from surrounding winds is also a possibility [6,18].
A turbulent fluid jet region can give rise to a variety of signals, from radio waves to γ rays. Furthermore, cascades of high-energy particles in the jets produce different particle populations linked by the transport phenomenon. The emergence of neutrinos that then leave the system can be detected by modern arrays.
In this work, the production of radio synchrotron emission and very-high-energy (VHE; up to hundreds of TeV) neutrinos from generic MQ jets is studied using the method of dynamic and radiative relativistic MHD simulations, where the model space is divided into computational cells. The effects of wind from a companion star are also included.
Adopting the one-zone (homogeneous) emission model in each eligible hydrocode cell (in this work, hydrocode refers to the PLUTO program [19]), the solution of successive transport equations connecting particle distributions in a cascade provides the emitted neutrino intensity as a function of local dynamic and radiative jet parameters. Particle cascade calculations, acting on hydrocode data outputs, employ the results from Monte Carlo simulations [20,21]. This way, a cell’s physical parameters are directly connected to its particle emission. The cascade timescale is assumed to be much smaller than the dynamic timescale of the modeled system. Repeating the above calculations over a number of different energies provides a neutrino energy spectrum at each jet point. Line-of-sight integration follows, leading to the production of synthetic neutrino images and model system spectra.
Charge neutrality is assumed in the jets, coupling the bulk flow proton and electron distributions dynamically through a relativistic MHD simulation. The latter provides the bulk proton density and also the magnetic field in each cell. High-energy hadron and lepton populations, obtained through shock front acceleration in the jets using the one-zone model in each computational cell, lead to synchrotron emission [6] taking into account local magnetic fields. Here, we focus on the radio synchrotron band, where the spectrum is flat or inverted and the emission from protons is essentially negligible (see Figure 5 in [6]). As shown in the Appendix A (Appendix A.1), the inner jet region of ≃10 3 pc is modeled here (corresponding to roughly 10 mas, at a distance to Earth of 5 kpc), where a flat radio spectrum typically occurs.
Lepton and hadron contributions that are comparable to radio synchrotron only occur when much more energy is assigned to hadrons than to leptons. As an approximation, only synchrotron emission from accelerated jet electrons is modeled here using the formalism in [22]. Self-absorption is also included.
Consequently, local MHD jet properties, as provided by the MHD model, allow the calculation of local radio synchrotron emission and absorption coefficients. Down the pipeline, this may lead to synthetic radio maps of the model system. When combined, the study of both synchrotron and neutrino emission from the jets may offer an improved understanding of high-energy processes in the system.
In the current work, the finite nature of the speed of light is taken into account when viewing the model jet, which constitutes an improvement over previous work [23]. Imaging a relativistically moving macroscopical object opens the window to a rather unexpected and even strange world of peculiarities. The basic mandates of Special Relativity, regarding length contraction and time dilation, constitute the mere beginning in the quest for comprehension of the actual appearance of a fast-moving object [24,25,26,27,28]. An observer will see the object view affected by a number of relativistic distortions [29,30,31].
Electromagnetic emission transformation from the jet frame of reference to the Earth frame requires performing the Lorentz/Poincaré transform [18,29,32,33]. Applying the latter transform for imaging purposes aims to reconstruct what the observer will actually see. Relevant to this point, [34] argues about the important difference between vision and measurements in Special Relativity, presenting that difference in a geometrical manner.
Radiation emitted from a jet is therefore subject to relativistic effects [18,35], including time dilation, relativistic aberration, and frequency shift, collectively leading to what is known as Doppler boosting and beaming [29,33,36]. Aberration causes the fast-moving object to actually appear rotated to a stationary observer [25,26,29,37], a phenomenon sometimes called the Terell–Penrose rotation.
Ray-tracing methods provide excellent-quality relativistic images, despite lacking in terms of efficiency compared to such techniques as polygon rendering [31]. In this work, a relativistic imaging method is employed, whereby time-resolved hydrocode data are crossed by rays called lines of sight (LOSs), either focused or parallel to each other. Each ray is then subject to beaming effects, which means Doppler boosting for E/M radiation, depending on the ray and the local velocity relative orientation.
Furthermore, the aforementioned high-energy proton distribution is Lorentz-transformed to a stationary frame, exhibiting dependence on speed and orientation, favoring particle emissions along the direction of the local velocity [38]. The resulting expression demonstrates relativistic boosting properties in a manner broadly comparable to the Doppler boosting of E/M radiation. A cascade of particle distributions then emerges, affected by local MHD properties, eventually leading to neutrino emission, which then escapes the system. Subsequently, with neutrinos, imaging “rays” move at the finite speed of light.
Simulating the above processes may help clarify the inner workings of the jets and their environs, leading to a more accurate description of the system of interest.
In the remainder of this paper, the theoretical background is presented for both radio synchrotron and neutrino particle emissions (Section 2). Section 3 briefly describes the employed software pipeline. The model setup is presented next (Section 4), including various model parameters. In Section 5, the results are presented and discussed. Finally, in Section 6, useful conclusions are drawn from the current work and possible future applications are proposed.

2. Theoretical Background

A conceptual link between electromagnetic and particle emissions from the jets may be formed. Proposed neutrino emission from jets favors the presence of high-energy protons and electrons, triggering particle cascades that lead to neutrinos. Jets are therefore acceleration sites, possibly at shock fronts inside them. High-energy electrons lay the foundation for synchrotron radio wave emission, as well as emission in other parts of the electromagnetic spectrum. Shock front acceleration therefore energizes the jets, leading to both particle and radiation production thereafter.

2.1. Radio Synchrotron Emission and Self-Absorption

Adopting high-energy electron and proton acceleration at shock fronts [22], a power-law high energy distribution is assumed for each of the aforementioned particle species. Electrons are discussed in this subsection; protons are discussed in the following one. For electrons,
N ( E ) = N 0 E γ c = K pl ρ 4 d E γ c ,
where K pl is a power law constant for high-energy jet electrons. K pl ρ 4 d = N 0 ; therefore, the reference electron density N 0 is taken as proportional to the local bulk flow proton density ρ 4 d based on charge neutrality. ρ 4 d is provided for spatiotemporal points by MHD code. Thus, a link is established between PLUTO’s hydrodynamic quantities and electron populations in the jet.
For accelerated protons and electrons, a high-energy cutoff is adopted according to the Hillas criterion [39]. E max = ZqB R s is the maximum energy E max achieved through acceleration for a particle of charge q and atomic number Z within a magnetic field B in an accelerator of size R s .
For the radio-scale model, the computational cell size of 10 13 cm is selected in order to fit the radio-resolved portion of the jet into the model space. Furthermore, the magnetic field is set at a value less than equipartition, at 10 G. The above equation will then result, for electrons, in a maximum energy on the TeV scale, well above the radio band which is of interest in this subsection. Furthermore, the energy cutoff from dominant cooling mechanisms [6] for radio-band synchrotron from a microquasar lies above radio wave energies. Thus, no cutoff is needed for electrons here.
On the other hand, for high-energy protons, a cutoff of E max = 10 6 GeV is employed [6,23].
Synchrotron emission and absorption coefficients are then calculated, following [22]:
ϵ ν = f pp ( 4 d ) K pl c 5 ( γ c ) ρ 4 d ( B 4 d sin ( LOS , B 4 d ^ ) ) ( γ c + 1 ) 2 f obs 2 D 4 D local c 1 ( 1 γ c ) 2 D 4 D local 2
k ν = f pp ( 4 d ) K pl c 6 ( γ c ) ρ 4 d ( B 4 d sin ( LOS , B 4 d ^ ) ) ( γ c + 2 ) 2 f obs 2 D 4 D local c 1 ( 4 γ c ) 2 D 4 D local 2 ,
where f pp ( 4 d ) is a proton profile function, with a threshold velocity u th of 0.1 c (f pp ( 4 d ) = 1.0 if | ( u / u th ) | > 1, and f pp ( 4 d ) = | ( u / u th ) | z fpp if | ( u / u th ) | < 1); f obs is the observing frequency; and f obs D 4 D local is the calculation frequency (blue-shifted or red-shifted according to the relative orientation between the local LOS and local velocity); B 4 d is the local magnetic field vector; LOS is the line of sight (LOS) vector at an angle ( LOS , B 4 d ^ ) to the magnetic field spatial vector; D 4 D local is the local Doppler factor of E/M emission from a computational cell, calculated as D = 1 u 2 ( 1 u cos ( losu ) ) , losu = ( LOS , u ^ ) being the angle between the LOS and the local velocity u, 0 ≤ u ≤ 1 (Appendix C.3.3); and γ c is the electron power law index, taken as 2 in the current work, equal—for simplicity—to the proton power law index. According to [22], the quantities c i , i = [1, 6], are c 1 = 6.27 × 10 18 , c 2 = 2.37 × 10 3 , c 3 = 1.87 × 10 23 , c 4 = 4.2 × 10 7
c 5 = 0.25 c 3 Γ 3 γ c 1 12 Γ 3 γ c + 7 12 γ c + 7 3 γ c + 1
c 6 = 1 32 3 × 10 10 c 1 2 c 3 γ c + 10 3 Γ 3 γ c + 2 12 Γ 3 γ c + 10 12 ,
where Γ is the Gamma function.
The above coefficients are computed at each point (computational cell) of the model jet system. The angle ( LOS , B 4 d ^ ) between the LOS and the local magnetic field B is calculated at each jet point, Appendix C.3.4.
In the above, only two Doppler factors and a frequency shift factor are included. The reason for this is to allow detection of the receding blob in the synthetic radio images. Additional Doppler factors can be employed, see Appendix C.3.

2.2. Nonthermal Proton Density

Neutrino emission from model jets comes from proton–proton interactions between a distribution of high-energy (non-thermal) protons and bulk flow protons [2,5,6,9,20,21,23]. Neutrino emission also comes from p γ interactions. For the energy range employed here, from 1 GeV to 5 × 10 5 GeV, the pp contribution is higher than the p γ one. According to Figure 7 in [6], at 1 GeV, the difference is several orders of magnitude in favor of pp, while even at 10 6 GeV, the pp contribution remains higher. Therefore, over the energy range of interest, the p γ contribution may, in the first approximation, be left out. Nevertheless, including the p γ contribution is a priority in future use of the model.
Some thermal protons gain energy at shock fronts within the jet according to the first-order Fermi acceleration mechanism within a time frame of [15,16,40]:
t acc 1 η c e B E p ,
where B is the magnetic field, E p is the non-thermal proton energy, e is the particle charge, and c is the speed of light. η is an acceleration efficiency parameter near the base of the jet for the process of efficient particle acceleration in moderately relativistic shocks [40]. In general, η depends on the shock velocity and on the local diffusion coefficient [41]. As an approximation, η = 0.1 [40].
A power law energy E distribution is employed for non-thermal protons, N p = K N p ( 0 ) E α [18]. The constant K provides the ratio between the number densities of hot and cold protons, and is much smaller than unity ( N p is the hot proton number density and N p ( 0 ) is the cold proton number density). In the model, the proton spectral index in the local jet matter frame, α ≈ 2 [5]. In general, α is a parameter of the model [9] that depends on the non-thermal proton distribution profile, and may be changed in future work.
In this work, transport equations are resolved for pions and neutrinos, while a power law is used for hot protons. Nevertheless, a transport equation may also be used in order to find the hot proton distribution [5].
The high-energy proton distribution is considered isotropic in the jet frame. The hypothetical anisotropy of the hot proton distribution can be reflected in the neutrino distribution [42], potentially introducing anisotropy in the jet frame’s particle emission field. Assuming that l sc < l r at each jet point, where l sc is the scattering length and l r the radiative length. The above approximation is based on the need to preserve, after every bounce, at least some proton energy [16]. The scattering length l sc is less than the radiative length l r , otherwise the proton would not have any energy left after the bounce, negating the acceleration process.

2.3. Neutrino Emissivity

For each computational cell, a high-energy proton distribution is transformed [38] from the cell’s frame to our frame, using the angle of local velocity to the LOS crossing that cell. In each cell, a local particle cascade emerges. From protons to pions and then to neutrinos, successive particle populations are linked by transport equations. As an approximation, the branch leading to muonic neutrinos is not included here. In each cell, transport equations are resolved along the particle cascade. Starting from a power law high-energy proton distribution, we obtain a pion distribution and then a neutrino population, which then leave the system [5,6,20,21]. The resulting neutrino population is
Q π ν ( E ) = E E max d E π t π 1 ( E π ) N π ( E π ) Θ ( 1 r π x ) E π ( 1 r π ) ,
where E is the neutrino energy, E π is the pion energy, N π is the pion number density, x = E / E π , r π = (m μ /m π ) 2 (m μ and m π are the muon and pion masses, respectively) and t π is the pion decay timescale. Θ ( χ ) is the theta function [6,43]. As shown in [5], E max = 10 6 GeV. The calculation leading to the above result can be found in [23] (see also the Appendix D). Neutrino emission is calculated for each eligible hydrocode cell. The imaging process may incorporate either parallel LOSs or a focused beam, where each LOS follows a slightly different path to a common focal point [44]. A synthetic image of the model system is thus produced.
Muon neutrinos are not included in the first approximation, since their anticipated contribution is perhaps an order of magnitude lower. Their inclusion is deferred to future work.

3. Computer Programs Used

3.1. RLOS

rlos [44] is an evolution of classical imaging code used in earlier works [45,46,47]. A ray, or an LOS, emanates from each pixel of the imaging side of the Cartesian 3D computational domain, Figure 1, or from a focal point, aiming at a point on a fiducial screen, Figure A4 (Appendix C). Either way, there is an imaging plane. Along the LOS, the radiative transfer equation is solved in each cell using local emission and absorption coefficients. Depending on the modeled situation, the coefficients may either be directly calculated or outsourced to another program.
rlos is organized into two outer spatial loops running over the imaging plane and an inner one-dimensional spatial loop, advancing in pairs of steps, one for each direction angle (azimuth and elevation), running over the length of an LOS (Figure A6). At the innermost point of the nested loop lies a conditional temporal loop, running over the hydro data time span (see also the Appendix C, Appendix C.6.1). Since the emission coefficients’ calculation load is global, it is performed, where feasible, before the loops in array-oriented operations in order to improve performance.
Lines-of-sight are drawn starting from a focal point (focused beam) or from a pixel of the yz or xz side (parallel rays) of the domain, Figure 1 (see also Figure A1 and Figure A4 in the Appendix C). Tracing their way along the given direction, they reach a length of ( x max 2 + y max 2 + z max 2 ) , where x max , y max , z max are the cell dimensions of the computational domain. In practice, on reaching the ends of the domain, an LOS calculation halts, and some LOSs may thereby end up being shorter than others. The above process is repeated within a 2D loop running over the imaging plane, with each LOS corresponding to a single pixel of the final synthetic image. As an approximation, no sideways scattering is considered along an LOS.
A model astrophysical system geometry may be directly inserted into rlos. As an example, a “conical” jet setup [48] is available to the user. Alternatively, data output from a hydrocode may be employed, as in the current paper, using PLUTO [19].
Figure 1. Three-dimensional schematic view of rlos applied to a model astrophysical jet for a parallel ray setup. The imaging side of the computational box is the yz plane located on the side of the box that appears closer to the reader. Lying on the yz plane, O is the point of origin of a random LOS with its own dashed coordinate system x y z . Alternatively, the imaging plane may also lie on the xz side of the box.
Figure 1. Three-dimensional schematic view of rlos applied to a model astrophysical jet for a parallel ray setup. The imaging side of the computational box is the yz plane located on the side of the box that appears closer to the reader. Lying on the yz plane, O is the point of origin of a random LOS with its own dashed coordinate system x y z . Alternatively, the imaging plane may also lie on the xz side of the box.
Galaxies 11 00110 g001

3.2. The PLUTO Hydrocode

PLUTO [19] is an open-source, 2D/3D modular hydrocode—a finite-volume/difference shock-capturing program—meant to integrate a set of (time-dependent) conservation laws. Initial and boundary conditions are conveniently assigned through an equivalent set of primitive variables. The relevant systems of equations may include hydrodynamics (HD), magnetohydrodynamics (MHD), and their special-relativistic counterparts, RHD and RMHD, respectively, in either two or three spatial dimensions. The solution of conservation laws is produced through discretization on a structured mesh—a logically rectangular grid surrounded by a boundary with additional ghost cells—in order to implement boundary conditions. The grid may either be static or adaptive, and various coordinate systems are available. The program may run efficiently in parallel on various platforms.

3.3. Nemiss

Nemiss [23,49] calculates neutrino emission and spectra from the output of PLUTO hydrocode. It stands between PLUTO and rlos, taking the burden of global particle cascade calculations off the shoulders of rlos, helping form the PLUTO–nemiss–rlos pipeline. Synthetic neutrino images are produced taking into account the finite speed of emitted neutrinos. Doppler boosting and frequency shifts are switched off in rlos when imaging with neutrinos.

3.4. Software Information

Intensity plots of the jet pair were created using Veusz (https://github.com/veusz), a software for plotting data written by Jeremy Sanders and contributors and distributed under the GNU/GPL license. rlos and nemiss [49], written by the author, are available under the lGPL license. PLUTO was written by Andrea Mignone and collaborators, and is available under GNU/GPL.

4. Model Setup

The MQ system is represented on two different scales: one for modeling neutrino emission and another for radio emission. On the smaller scale, meant for neutrinos, an accretion disk is assumed around the collapsed object [50], while the companion star itself lies outside the model space. A continuous ejection, representing the beginning of a new blob, is employed.
On the other hand, on the radio emission scale, both participants of the binary system lie within the same computational cell, while a sequence of plasmoids is employed.
In both cases, twin jets emerge from near the compacted star, with a collimating toroidal magnetic field component. The field is set higher at the neutrino emission scale.

4.1. Radio Synchrotron Emission Model

A twin model jet system is synthetically observed using radio synchrotron emission at a “typical” radio frequency of 8 GHz. A 3D homogeneous Cartesian coordinate system is employed. Jets enter the grid at a generic speed of 0.8c (a speed attributed, for example, to the jets of GRS1915+105), emanating from a central location in the grid.
The jet is synthetically observed using synchrotron emission and self-absorption [22]. In each computational cell, the angle between the local magnetic field and the LOS is calculated, leading to determination of the cell’s emission and absorption coefficients. Likewise, the local angle between the LOS and velocity enables calculation of each cell’s Doppler boosting.
A sequence of twin relativistic blobs is ejected from the jet base, moving in opposite directions at an angle to the observer. One jet is approaching, the other is receding. The model jet is made of a series of such plasmoids ejected for 30 hydrocode time units every 100 time units. The latter choice of the plasmoid size is aimed to produce typical intermittent jets made from a series of blobs. Two separate hydrocode runs were performed, one with a heavier jet and another with a lighter jet.
A special note should be made concerning the employed model jet power. The heavier radio jet in the model is very strong, at about 10 43 ergs 1 , whereas the lighter one stands at around 10 41 ergs 1 . These values are higher than the Eddington luminosity for an MQ system, especially the heavier jet model. The reason for this is that limited computational resources and the use of a homogeneous grid led to a low model resolution being employed. The latter defines the minimum nozzle diameter, which is much larger than the real jet. Thus, in order to keep jet densities realistic, the jet powers tend to be higher than normal. Details of the runs are shown in Table 1.
A toroidal magnetic field of 10 G is introduced, resulting in a magnetic energy density below the kinetic energy density (see the Appendix B). For simplicity, the magnetic field in this run was initially wrapped around the y-axis.
At regular time intervals, a model system instance was saved to disk. Then, synchrotron emission ϵ ν and absorption κ ν coefficients are calculated at each spatial point for a given set of observation angles within rlos. A series of four-dimensional arrays is thus obtained, through which lines-of-sight travel. Rays start at the moment of observation and move backwards in both space and time in order to meet blobs at an earlier instant.
The angle to the y axis was 55 degrees in order to enhance the effect of the apparent acceleration of approaching plasmoids in the synthetic radio images. The finite ray speed focused beam imaging method was employed in all radio imaging runs. As a test, an RMHD run was also performed with an angle to the LOS of 55 degrees in order to better visualize apparent superluminal motion. A simplified stellar wind construct was employed; its density was inversely proportional to the distance from the companion star [45].

4.2. Neutrino Emission

In the neutrino-scale model scenario, twin hadronic jets were simulated with PLUTO code [23]. Jets are viewed from the side, while the finite-ray speed imaging mode of rlos was employed.
For the given model view orientation, neutrino emission was calculated at each eligible spacetime point of the computational grid using nemiss [49]. The neutrino emissivity is separately calculated in individual computational cells using the angle (los,u) formed between the LOS and the local velocity. The parameters used in this scenario are in Table 1. The binary companion now resides outside the computational grid at the position (400, 0, 400), a binary separation distance of 4 × 10 12 cm, affecting the model with its stellar wind [23], which falls off as 1/r 2 away from the star, a typical setup for stellar wind. The accretion disk was approximately simulated as disk-shaped, and a basic accretion disk wind construct, falling off as 1/y away from the disk, was also included.
The jet being continuous is a feature compatible with the radio-scale model, since at the ν scale, the time unit of the model is 1000 times less than in the radio scale. The beginning of the injection process of a single radio plasmoid was conceptually modeled in the neutrino-scale scenario.
Furthermore, rlos [44] was employed, which reads the combined results of PLUTO and nemiss, producing synthetic neutrino images of the model system using the focused beam geometry setup.
At the ν scale, the relativistic transformation of the hot proton distribution in [38] was employed. The magnetic field is toroidal and was adjusted for equipartition, B r z = 8 π ρ r z [15,16] (Appendix B). The above choice of magnetic field is an approximation suitable for the inner jet region, with strong fields threading the jets. External magnetic fields can be taken as smaller [51]; therefore, as a first-order approximation, they are omitted from the model’s surrounding winds.
The neutrino emission in each computational cell was calculated using the formalism presented earlier in this paper (Section 2.3). This method, albeit costly, allows one to obtain separate neutrino emission from each spatiotemporal point of the model. Thus, we aim for the result
I ν = I ν ( r , t )
where the 3D space location r is represented by the x, y, and z coordinates of the computational cell in question. Time t is obtained from the time tag of the hydrocode snapshot to which the cell belongs (MHD datasets are four-dimensional, including the dimension of time). The above equation is globally applied to all selected PLUTO data (the user may select beginning and end times for the global calculation).
For reasons of economy, a double filter was applied, whereby neutrino emission was calculated only in cells in which the velocity is not too far from the LOS (cos(losu) greater than 0.08) and whose speed is at least 0.1 c.
On the neutrino scale, the resulting jet’s kinetic power was L k = 2 × 10 38 (see the Appendix A). In [6], the authors argue a 10% Eddington jet kinetic luminosity, leading to L k = 10 38 ergs 1 for a 10 M black hole, which is comparable to the current case. In [6], either L p L e 100 or ≃1 was adopted; we selected the former, which implies a hadronic jet. As an approximation, neutrino emission from a high-energy proton distribution is obtained, omitting effects (such as synchrotron emission cooling) from the corresponding high-energy electron distribution in the jet. The jet base is situated near the center of a Cartesian grid. The same spatial resolution is employed on both the radio and neutrino scales, resulting in a higher jet power in the radio-scale models. Nevertheless, this can be offset by a normalization process, while keeping densities, which strongly affect the emission and absorption calculations, realistic.
The current work constitutes an improvement over previous work [23]. More specifically, time synthetic neutrino images are produced, leading to a more detailed view of the model system, as opposed to merely synthetic energy spectra in [23]. Furthermore, the finite speeds of neutrino “rays” are now incorporated in the model, leading to more realistic neutrino imaging of model jets. What is more, the current work focuses on sidereal emission from model jets, whereas [23] provided a more general study of neutrino emission from jets viewed at different angles.

4.3. Model Parameters

Table 1 shows a number of simulation parameters. These include the computational cell length, jet density, and wind maximal densities (which gradually decline away from their sources).
On the neutrino scale, the jet kinetic luminosity is 2.5 × 10 38 erg/s, using a low spatial resolution of 60 × 100 × 50. Due to limited computing resources, such a low resolution was necessary in order to accommodate both the heavier neutrino emission calculation that followed and the time-resolved nature of the calculations, which essentially leads to four-dimensional datasets.
In the radio scale model, the use of temporally resolved data, consuming computing resources, and the use of a fixed homogeneous grid, also necessitates a smaller spatial resolution, namely 60 × 100 × 50, and thus a larger computational cell. Radio-scale model jets are then more powerful, in order to keep their densities realistic, at the low resolution employed. At the jet nozzle, this means a thicker jet. Opting to keep jet density realistic, the overall kinetic luminosity then becomes higher than normal, mainly in the heavy jet radio model. Nevertheless, a normalisation process (electron power law constant K) may, to a certain extent, absorb the effects of an artificially increased jet power. On the other hand, synchrotron emission and absorption coefficients do depend, among other things, on the local density, so a need arises to focus on the density in this context.
In PLUTO, the piecewise linear method was set up using the MUSCL Hanckock integrator with an ideal equation of state. In the neutrino-scale simulation, the binary companion is located outside the grid, and was estimated to be at most up to an order of magnitude larger than the compact object. As mentioned above, the jet speed is 0.8c.

5. Results and Discussion

5.1. General

The twin jet simulations in this paper represent three fiducial microquasars, dynamically set up to resemble a number of actual microquasars. In the synthetic imaging process (rlos code), a focused beam geometry was employed in combination with a finite ray speed. Synthetic images were projected onto a fiducial imaging screen, parallel to the side of the grid (YZ), before the beam converged to a focal point, Figure A4.
In general, the imaging process may or may not use all snapshots available to it depending on the light crossing time of the model segment (adjusted through the clight parameter in rlos, Appendix C.4). Trying to read more snapshots than are loaded corrupts the hydrocode time array of rlos, called T, resulting in errors.
For the neutrino calculation, as mentioned above, a double filter was used for the velocity and for the (los, u) angle. A minimal velocity and a maximal angle were set in order to trigger the neutrino emission calculation for a particular cell. This way, the expensive part of the simulation was only performed when it was really worth it. This partly alleviated the discrepancy between computational costs of the dynamic and the radiative parts of the model.
For the radio calculation, a similar filter was applied for the (los, B) angle. The use of filters highlights sideways emission from jet elements, with near-relativistic velocities roughly perpendicular to the jet axis.
An important aspect of this modeling approach is that a cell has a different visible emissivity from Earth to that of its neighbors. This is because each cell may differ from the next one in terms of velocity value and orientation to us. Consequently, the output of a hydrodynamic model will differ from that of a steady-state model. Individual flow elements may appear either boosted or de-boosted depending on their velocity’s orientation to the observer.
Even from the side, adequate emission may be obtained due to elements of the flow moving relativistically sideways, especially during the initial jet front expansion. Furthermore, delays in the arrival of emission emanating from the inner part of the jet mean that the initial violent interaction between the jet and the surrounding wind still affects synthetic images taken at subsequent time instants. Inner winds that are heavier than the jets, as was the case in the model runs, further prolong the effects of sidereal emission on the images by temporarily constraining the jet’s “bubble” near the jet base. Thus, initial jet expansion still echoes in images taken later on in the simulation.
The distances in the synthetic images are in computational cells, where, in this work, one cell = two hydrocode length units. On the other hand, in the hydrocode data renderings, the scale is in hydrocode length units. For each scenario, the hydrocode length unit is shown in Table 1.

5.2. Radio-Scale Model Results

In the radio-scale model, ejected plasmoids form an intermittent twin jet (Figure 2 and Figure 3) and inflate a twin cocoon while traversing the companion star’s stellar wind.
Moving away from the jet base, plasmoids seem to expand while crossing ambient wind, thanks in part to a declining wind density away from the binary system. A relatively low toroidal magnetic field further facilitates plasmoid expansion in the model.
An equatorial concentration of wind, and possibly jet, matter forms dynamically (Figure 2 and Figure 3), leading to some emission from that region, as seen in the corresponding synthetic radio images that follow.
We can observe the apparent acceleration of the approaching plasmoid (Figure 4 and Figure 5) on the sky plane (fiducial model screen), while the receding blob moves slower. The approaching blob is also brighter than the receding one. What is more, taking into consideration the finite nature of the speed of light, earlier dynamics affect images taken at later times. Synthetic radio images of the model system in general exhibit a delay versus actual hydrocode plots bearing the same time tag (Figure 4 and Figure 5).
A toroidal magnetic field threading the blobs leads to synchrotron emission in the direction of the observer, even from the receding blob (only two Doppler factors are employed to allow some visibility of receding plasmoids). The model ejected plasmoid pair exhibits a quick rise in intensity and a more gradual decrease over time, Figure 6. In Figure 6, for the heavy jet case (top curve), a previous pair of blobs happens to leave the computational grid at the time of new blob pair insertion, so a drop in intensity occurs there. This drop is observed at model time 150, but it actually took place shortly before the end of injection of the pair of blobs (since there is a recovery in intensity after the drop and thus injection continues after the drop), at around time 120–125 (injection of a second blob pair ends at model time 130).
A suitable jet orientation setup has the potential for apparent superluminal motion of the approaching jet. This is visible in Figure 4 and is particularly notable in Figure 5.
As a test, a simulation was also run with an approaching “rectangular” blob. A frontal synthetic image of the blob, moving straight into the fiducial observer (Figure 7), demonstrates the effects of a finite ray speed, as central rays from the rectangle appear stronger than those from its corners due to the delayed arrival of the latter.

5.3. Neutrino-Scale Results

The PLUTO hydrocode was run in order to simulate the jets on the neutrino-emission scale (Figure 8). A number of auxiliary PLUTO user parameters were also defined beforehand in order to accommodate particle emission results from each cell later on. The purpose was to prompt PLUTO to create additional data files filled with a random number. These additional data files, while filled with meaningless data for now, still form part of the PLUTO data save system. They are meant to act as a vehicle in order to accommodate particle emission results later on, with one file for each neutrino energy.
The nemiss program was then run to calculate neutrino emissions from hydrocode data for a specific imaging geometry and setup. This program is able to read 4D spatiotemporal data outputs from PLUTO and convert them into a 5D array, including the particle energy as the fifth dimension. Then, nemiss calculated the neutrino emission at each point of the 5D data array. Nemiss results were then saved into the aforementioned suitably prepared data files of the PLUTO hydrocode. Thus, nemiss processes the PLUTO output to add a neutrino emission spectrum at each spatiotemporal data point.
PLUTO data processed by nemiss were then ready to be read by the relativistic imaging code rlos [44], which produced synthetic neutrino images of the system. The sum of intensities over all synthetic jet images was calculated over a range of particle energies, leading to a plot of jet neutrino intensities, Figure 9. Furthermore, the model was run at two different time instants: shot number 45 (t = 90 s) and shot number 90 (t = 180 s). A relativistic imaging process was used to produce synthetic images, Figure 10.
Figure 8 shows narrow jets that are slowly expanding into the surrounding winds, collimated by a strong toroidal magnetic field component (Figure 11). This small half-angle is then rather counter-intuitively expected to result in a faster decline in neutrino emission with energy, as discussed in [6]. Jets interact with winds, their cocoon expanding sideways as well as forwards into the surrounding wind matter. Furthermore, sidereal expansion of the inner part of the cocoon accrues enough dense and relatively fast matter to achieve adequate (beamed, according to [38]) sideways neutrino emission, as seen in the synthetic neutrino images that follow. The effects of the accretion disk and its wind construct therefore seem to play an important dynamical role in providing sideways particle emission from the jet system. This is important for sideways neutrino emission (Figure 10) because it shows localized neutrino emission along a direction perpendicular to the jet axis.
As mentioned earlier, adequate neutrino emission may occur towards Earth, even from MQ jets not aligned with the LOS to Earth. This leads to a rather increased number of MQ candidates for neutrino detection, particularly those with rich dynamic interactions with the surrounding winds.
The above result positively affects emission from a galactic MQ distribution [23]. Even in MQs whose jets point perpendicularly to the LOS to Earth, some relativistic boosting may still appear in parts of the flow moving towards us, especially early on in the ejection event.
As shown in [23], the scale of total emission is expected to increase the closer the LOS is to the approaching jet axis. The current results (Figure 12) demonstrate the possibility of potential observations, even when the jet is observed from the side. In other words, in this work, beaming [38] is localized in cells whose relativistic speed points towards Earth, even though the jet axis lies essentially perpendicular to our line-of-sight. As a concrete example, the detection ability of a cubic km array is depicted in a normalized SED plot, Figure 12, presenting a reasonable possibility for detection.
The above estimate may then be employed to provide a rough estimate of expected neutrino emission from a microquasar distribution in the galaxy. The authors of [52] proposed an estimated population of around a hundred systems in our galaxy. Furthermore, their discussion of γ ray emission from microquasars clarified the importance of relativistic boosting in jet emission. Thus, orientation to Earth plays a major role here, and the situation is similar for neutrino emission [38].
For neutrino emission, in a manner similar to [23], we proceed by assuming 100 systems at various distances ranging from a minimum of 1 kpc to a maximum of 30 kpc, with an average kinetic luminosity similar to our model system. The linear dependence of emission on the latter quantity facilitates such a simplification. A distance of 1 kpc commands a flux at Earth of 25 times more than our model value, whereas the representative system situated at 30 kpc has 36 times less than at 5 kpc.
In this analysis, the angle of the line-of-sight to the x axis is employed to ensure compatibility with the rlos geometry setup. This angle is complementary to the angle to the jet axis.
An orientation of less than 60 degrees might lead to an intensity an order of magnitude lower than our value, but a jet system aimed towards us could perhaps have up to 100 times more visibility on Earth [23]. The latter point can be somehow altered by the present analysis insofar as, even from the side, the expanding cocoons lead to local sidereal beaming towards Earth.
Figure 12. Neutrino-normalized intensity spectral emission distribution plots at snapshots 45 and 90 (t = 90 s and 180 s, respectively). The synthetic imaging viewing angle of the model system used to produce these spectra is from the side, roughly perpendicular to the jets. The contribution to the emission is mostly from matter elements moving sideways during jet cocoon expansion through ambient winds. Sizeable sideways emission does appear in the synthetic images, and that is reflected in the normalization process assumptions used here. It should be noted that the quantitative intensity difference between the earlier and later snapshots is due to the actual dynamical evolution of the jet system in the PLUTO RMHD simulation, as both spectra were produced using the exact same normalisation assumptions. The two spectra exhibit similarities, as expected by the similarity of the synthetic images in Figure 10. We can also see the reference sensitivity of a cubic km detector array in the first approximation (dotted line), and an improved version depending on the energy (thick line) [53].
Figure 12. Neutrino-normalized intensity spectral emission distribution plots at snapshots 45 and 90 (t = 90 s and 180 s, respectively). The synthetic imaging viewing angle of the model system used to produce these spectra is from the side, roughly perpendicular to the jets. The contribution to the emission is mostly from matter elements moving sideways during jet cocoon expansion through ambient winds. Sizeable sideways emission does appear in the synthetic images, and that is reflected in the normalization process assumptions used here. It should be noted that the quantitative intensity difference between the earlier and later snapshots is due to the actual dynamical evolution of the jet system in the PLUTO RMHD simulation, as both spectra were produced using the exact same normalisation assumptions. The two spectra exhibit similarities, as expected by the similarity of the synthetic images in Figure 10. We can also see the reference sensitivity of a cubic km detector array in the first approximation (dotted line), and an improved version depending on the energy (thick line) [53].
Galaxies 11 00110 g012
Orientation then remains the most important factor, but the low intensity end of the range might be updated upwards, since the local cell velocity might still point towards Earth even if the jet axis does not. Then, we have the distance, and lastly, the jet kinetic power.
The above order of importance allows for an estimate of perhaps 5% or, in other words, five systems with a very high relativistic boosting towards us. Then, we have maybe 40 or 50 sources aimed at angles above 45 degrees but below 90 degrees, and lastly maybe 50 sources at below 45 degrees which still may have quite sizeable contributions. The high boost sources probably contribute the most on average, and the ones viewed from the side have a smaller effect. A possible system located at a smaller distance, oriented towards Earth, would of course dominate the above distribution, but the possibility of such an occurrence is not very large.
On the basis of the above discussion, we then use a rough average for a neutrino-emitting galactic microquasar located at 15 kpc, with the kinetic luminosity of our model (less affecting factor) and orientated at 37.5 degrees from the line-of-sight.
The reason for using an average angle of 37.5 degrees, higher than the estimate of 30 degrees in [23], is the higher contribution from systems not aimed towards us due to increased sidereal emission from dynamical jet systems.
The orientation of the local velocity and the magnetic field seems to play a crucial role in both neutrino and radio emission calculations. Consequently, the differential projection in each cell strongly affects the final images.
It seems possible, even more than in [23], that the detection of background emission from a potential distribution of microquasars in the galaxy lies within the realm of modern detector arrays. This is also a consideration for the next generation of new or upgraded arrays. On the other hand, a single X-ray binary system might act as a galactic source of high-energy neutrinos. This is a potential target for a particle sensor with an increased angular accuracy. The variability of microquasars within the human timescale, combined with their relative stability as a known point source, offers a good target for observation, especially combined with sensors operating in the electromagnetic spectrum.

6. Conclusions

In the radio-scale model, apparent acceleration of the approaching plasmoid is observed on the fiducial sky plane, as well as an increased brightness. On the other hand, the receding plasmoid appears to move slower, while also appearing dimmer. In the RMHD hydrocode model, each blob pair is essentially identical, so the imaging model can capture the apparent relativistic acceleration and beaming at each jet point. Furthermore, the frequency shift is also included separately for each computational cell. This way, a more realistic simulation of the model system is achieved, facilitating a comparison with observations. The dynamics in the hydrocode and the magnetic fields threading the jets are then better connected to the actual appearance of the system to a detector array.
In the neutrino-scale model, an updated result is obtained, favoring sideways emission from relativistic jets. This may lead to an increased number of MQ candidate neutrino sources.
The synthetic imaging method using special relativistic methodology may be expanded to a general relativistic framework. At every grid point, a gravitational potential may alter the course of a ray, depending on the local matter and energy properties, as provided by a suitable hydrocode simulation. This way, a broader range of astrophysical problems can be approached with increased realism. Particle emissions can be included, if suitable transformation equations are provided for their energy spectra from the source reference system to the stationary system.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://github.com/teoxxx/nemiss_proper.

Acknowledgments

We are grateful to our colleagues for their valuable comments on the manuscript. Special thanks go to G. E. Romero (UNLP, IAR) for their suggestions on improving the content of this work.

Conflicts of Interest

The author declare no conflict of interest.

Appendix A. Jet Size and Jet Energetics

Appendix A.1. Model Space Size

The computational grid size is l grid = 200 × 10 13 cm = 2 × 10 15 cm. Distance to the generic MQ model system is taken as D = 5 kpc, or D ≃ 1.5 × 10 22 cm. Therefore, sin( l grid D ) ≃ 1.3 × 10 7 . For such small angles, sin(angle) ≃ angle. Thus, arc size is around 1.3 × 10 7 rad, or roughly 27 mas. The latter size is for the full computational domain, including jet and counterjet, as well as certain empty space margin. Consequently, for the core of the model grid, the model jet and counter-jet region spans no more than perhaps 15 mas. We may then consider a synchrotron emission region of up to 10 mas for a model jet, which is generally thought to be the inner compact region, where a flat to inverted spectrum may occur.

Appendix A.2. Normalization

The jet kinetic energy at its base is [6]
E k = ( Γ 1 ) m c 2
A constant jet speed u is considered for this calculation, meaning a constant Γ Lorentz factor. m is the mass of a jet portion traversing a jet cross section near the base. Then, jet kinetic power P k is the kinetic energy traversing the cross section per unit time
P k = d E k / d t = ( Γ 1 ) ( d m / d t ) c 2
where the speed is taken to be constant during an ejection episode (it was also set to be constant in the simulation described here). However,
d m / d t = ρ d V / d t = ρ A d x / d t = ρ A u
where ρ is jet density, V is the volume of a thin jet cross section, of length dx, near the jet base, and A is the jet base cross section area, also taken as a constant both in the simulation and here. Therefore,
P k = d E k / d t = ( ( Γ 1 ) ρ A ) c 2 u
or
P k = d E k / d t = ( ( Γ 1 ) ρ N cell L cell 2 ) c 2 u
where A = N cell L cell 2 is the area of jet cross section near the base, L cell is the length of the edges of a cubical computational cell, at the jet base, L cell 2 is the frontal area of a cell, and N cell is the number of such cells forming the jet base cross section.
We then express density as a function of proton number density N p and proton mass m p
ρ = N p m p .
The field of view (FOV) of the synthetic image is taken at an estimated 1/10 of the total solid angle. In-there, a similar portion of total emission is assumed, allowing for the majority of emission to go in the vicinity of the jet axis orientation. In total, a field-of-view factor of 1/100 of total emmission is employed: f fov = 0.01.
For the neutrino-scale model, let us define neutrino luminosity L ν as the power emitted through neutrinos from the jet, which is a fraction α of the total kinetic jet power (jet kinetic luminosity P k = L k ).
Consequently, α = L ν / L k , represents the portion of jet power (total power of hot and cold protons) emitted in neutrinos (the latter portion, α , is different than the portion of the hot proton energy converted to neutrinos, which usually bears a value of the order of a few percent). For normalisation, a working value is taken as α = 10 3 . This is compatible with a value of q rel = 0.1, for the energy portion, q, of hot protons in the jet [6,9]. It is out of energy portion q that a few percent end up to neutrinos. In summary, as an order of magnitude approximation, if jet kinetic energy is 100, then hot proton energy is taken as 10, and neutrino energy is taken as 1 percent of that 10, which is 0.1.
The shape of the spectrum is also affected by acceleration efficiency [9], and from the opening angle of the jet [6], thus affecting the area under the neutrino spectrum plot. As an approximation for the above effects, we adopted a value of 0.01 for the energy transfer from nonthermal protons to the neutrinos.
We also set u = β c. A less-than-unity positive filtering factor f f is employed that accounts for not using all jet cells, but only those with velocity orientation closer to the LOS and with speed above a given limit. We then have
L ν = α L k = α P k = α d E k / d t = f fov f f α ( Γ 1 ) ( N p m p N cell L cell 2 ) β c 3
The intensity of the jet is then expressed as I ν = L ν / 4 π D 2 , where D is the distance to Earth. Thus,
I ν = f fov f f 1 4 π D 2 α ( Γ 1 ) ( N p m p N cell L cell 2 ) β c 3
In the neutrino-scale simulation, the jet beam travels at β = u c = 0.8, with a density of 10 11 protons/cm 3 . L cell is 10 10 cm, while the number of cells comprising the beam at its base at this resolution is N cell 15. Distance to Earth is taken here with a typical value of D = 5 kpc or approximately 2 × 10 22 cm.
We then integrate the area under the curve of the arbitrary units plot (Figure 9) for our case of viewing the jet from the side. That case is supposed, for the purposes of normalization, to be the one matching the orientation of the hypothetical system in relation to Earth. We perform a cumulative sum over the roughly 10 points. Thus, we find about 10 11 , in arbitrary units (AU)*GeV. We replace an AU with a constant C 0 , so that AU = C 0 erg/(s*cm 2 ). We set I ν = L ν / 4 π D 2 equal to the area under the plot of Figure 9, called PLOTAREA, expressed in units of C 0 , in order to find the latter (normalization constant)
I ν = f fov f f 1 4 π D 2 α ( ( Γ 1 ) ( N p m p ) N cell L cell 2 ) β c 3 = ( PLOTAREA ) C 0 erg / ( s cm 2 ) GeV
For our case, we find C 0 ≃ 10 20 , which is the value of the arbitrary unit C 0 . Using the above constant, we multiply by it the value given in arbitrary units for the particle emission. Thus, the intensity plot is multiplied, and we arrive to the updated plot in Figure 12, which may be directly compared to other models and to observations.
For the radio-scale model, we may consider an estimate of the intensity in radio at 8 GHz. Relativistic beaming is now present, for the E/M radiation. In a manner similar to particles (see above), beaming leads to an approximate value for the field of view factor of f fov = 0.01. The portion of the jet kinetic power L k , emitted at the radio band in 8 GHz is estimated as PORTION.
For L k = 10 38 ergs 1 , L k × PORTION = (Power) 8 Ghz = 4 π D 2 × I ( Earth ( 8 GHz ) ) . Adopting 1 Jy (1 Jy = 10 23 ergs 1 cm 2 Hz 1 ) as broadband intensity I at Earth, and a channel width, at 8 GHz, of 10 2 GHz (10 2 GHz = 10 MHz = 10 7 Hz), I ( Earth ( 8 GHz ) ) = 10 23 + 7 = 10 16 ergs 1 cm 2 . Then, a channel power value of
( Power ) 8 Ghz = L k × PORTION = 10 38 ergs 1 × PORTION = 4 π D 2 × I ( Earth ( 8 GHz ) ) = 4 π D 2 × 10 16 ergcm 2 s 1 = 4 π × ( 4 × 10 44 cm 2 ) × 10 16 ergcm 2 s 1 10 30 ergs 1
is obtained. Consequently, the value of PORTION can be estimated to be of the order of 10 8 . Therefore, our image includes (1/30) × (PORTION) × (L k ).

Appendix B. Equipartition

The equipartition calculation follows. As shown above, jet kinetic power is
L k = P k = d E k / d t = ( Γ 1 ) ( N p m p N cell L cell 2 ) β c 3
where d m d t = ρ d V d t = ρ A d x d t = ρ A u
Kinetic energy density is [6]
ρ k = L k π R j 2 u j = L k A u = ( Γ 1 ) ( N p m p N cell L cell 2 ) β c 3 A u ,
At the nozzle, for the ν -scale model, ρ k 1 × 10 8 ergcm 3 . For the heavier jet radio-scale model, ρ k 1 × 10 9 ergcm 3 (3 × 10 8 ergcm 3 temporal average due to intermittent jet), and for the lighter jet radio-scale model, ρ k 1 × 10 7 ergcm 3 (3 × 10 6 ergcm 3 average).
Also
B = 8 π ρ B
and
ρ B = B 2 / 8 π
For both radio models, ρ B 10. For the ν -scale model, ρ B 10 7 . For equipartition, the kinetic and magnetic energy densities have to be equal to each other, ρ k = ρ B . Thus, both radio models have kinetic energy density higher than equipartition, whereas the ν -scale model roughly fullfills the equipartition assumption.

Appendix C. RLOS Special Relativistic Imaging Code

Additional properties of rlos imaging code are described here.

Appendix C.1. Theoretical Background for Imaging Code

Refs. [54,55] provide an early computerized attempt to reconstruct a relativistic image, through the eyes of an observer crossing a scene at high velocity. Ref. [56] demonstrates the importance of the relativistic transform of brightness and color. When imaging a jet, these correspond to Doppler boosting and frequency shift, respectively. [56] discusses an object that moves at uniform speed across the field of view, but is visually large enough for the angle between velocity and line-of-sight to vary along the object. Applying the Lorentz transform changes brightness and color in a separate manner, for each point of the observed object. Ref. [29] improves on such calculations, providing various methods for relativistic visualization, in both Special and General relativistic frameworks.
Ref. [57] calculate the visual appearance of wireframe relativistic objects, by mathematically inverting the course of light, from an image point to the emission event. They provide expressions that directly describe how a series of objects would look like, when moving at high speed, in front of a stationary observer. The efficiency of their method is then compared to the increased detail of a related ray-tracing project [58]. Ref. [36] image scenes with a fast observer traveling through their artificial environment. They also relate their simulations to actual imaging experiments, using the femto-photography technique [59]. Furthermore, they introduce a number of additional details into their models, such as camera distortions from traveling at very high speed. Ref. [60] present a framework, where the subject of relativistic imaging is explored, in an accessible manner.

Appendix C.2. Time-Resolved Imaging

Appendix C.2.1. Accessing 4-Dimensional Data

The finite nature of the speed of light crucially affects the appearance of a fast-moving object. Consequently, drawing a relativistic image of an astrophysical system, necessitates the availability of information regarding both its spatial properties and its temporal evolution. In the present case, when executing the hydrocode, before running rlos, the temporal density of snapshots to be saved to disk at regular intervals, was suitably adjusted. The smaller those intervals are, the better the temporal resolution of hydrocode data are. A series of snapshots are then loaded to RAM by rlos, which thereby requires more memory in order to run properly than the hydrocode itself does. Time is measured in simulation time units, which are read by PLUTO’s attached pload.pro routine, which loads data into rlos.
The total time span available to an LOS (as measured in simulation time units, not merely in number of snapshots), Δ t LOS ( total ) = t ( last shot ) t ( first shot ) (not to be confused with interval Δ t shot between successive snapshots) should be preset to be larger than light crossing time of the model system, for the selected LOS angle settings. Documenting model jet evolution generally requires hydrocode data saves to be rather dense in time, especially for fast-changing flows. On the other hand, a lower temporal resolution probably suffices for a steadier, slower-paced flow.

Appendix C.2.2. Traversing 4D Arrays

Introduction

A series of hydrocode snapshots are loaded to RAM, populating the elements of 4-dimensional (4D) arrays. Let us examine the case of backwards in time imaging calculation. From a temporal point of view, we begin from the simulation time corresponding to the last of the loaded snapshots, called shotmax. From a spatial point of view, we start at a point of the imaging plane, which is a side of the computational box (Figure 1). As the calculation advances, in 3D space along the LOS being drawn (Figure A1), the algorithm keeps checking whether to jump to a previous temporal slice while staying on target in 3D (Figure A2). Consequently, the LOS advances back in time through data by accessing different instants from the 4D data arrays. As a test case, the LOS may also be set to advance forward in time, beginning from the time tag of shotmin, the first of the loaded snapshots (Figure A2 and Figure A3).
Figure A1. Schematic of the spatial propagation of a line-of-sight (LOS) through a 3D Cartesian computational grid. In the discrete grid, according to the design of the algorithm, there are 3 available directions to be taken at each step along the LOS: right, up, and climb. These correspond to x, y, and z, respectively. During propagation, the LOS tries to follow its given direction, as defined by the two angles of azimuth and elevation. More specifically, every two steps, a decision is first made on azimuth, either right or up. Then, for elevation, it is either climb or another azimuth decision. Along the LOS, horizontal steps point to the ‘right’ direction. Diagonal steps represent going ‘up’, while vertical ones constitute ‘climb’ steps.
Figure A1. Schematic of the spatial propagation of a line-of-sight (LOS) through a 3D Cartesian computational grid. In the discrete grid, according to the design of the algorithm, there are 3 available directions to be taken at each step along the LOS: right, up, and climb. These correspond to x, y, and z, respectively. During propagation, the LOS tries to follow its given direction, as defined by the two angles of azimuth and elevation. More specifically, every two steps, a decision is first made on azimuth, either right or up. Then, for elevation, it is either climb or another azimuth decision. Along the LOS, horizontal steps point to the ‘right’ direction. Diagonal steps represent going ‘up’, while vertical ones constitute ‘climb’ steps.
Galaxies 11 00110 g0a1

Time-Resolved Imaging Calculations

For every LOS, there is a point of origin (POO), located on the imaging side of the computational grid (Figure 1). That point, addressed in rlos code as (nx10, ny10, nz10) and here as O , is the beginning of the LOS’s axes x , y , z , parallel to x, y, and z, respectively. A 2D loop covers the imaging surface, the POO successively locating itself at each of the latter’s points.
As we progress along an LOS, a record is kept of where we are, in 3D space. This record comprises the LOS’s own integer coordinates, rc, uc, and cc, measured in cells, starting from its POO. The above symbols stand for right-current, up-current, and climb-current, representing the current LOS advance in the x , y , and z axes, respectively (Figure 1 and Figure A1). The current ray position is then (nx10 + rc, ny10 + uc, nz10 + cc).
A timer variable, curtime (standing for current LOS time), is introduced for each LOS, recording the duration of insofar ray travel along the LOS. The aforementioned timer is preset at the beginning of each LOS to the hydrocode time of the first loaded data snapshot (forward in time integration), or of the last snapshot (back in time integration). For backwards in time ray-tracing, the above duration is subtracted each time from t(shotmax).
Figure A2. Three successive instants of a line-of-sight traversing a jet. At regular intervals, we jump to a new 3D slice of a 4D spacetime array, obtaining a discrete approximation of the time continuum in the form of hydrocode snapshots. The calculation may proceed either forward in time (from bottom to top) or backwards in time (from top to bottom).
Figure A2. Three successive instants of a line-of-sight traversing a jet. At regular intervals, we jump to a new 3D slice of a 4D spacetime array, obtaining a discrete approximation of the time continuum in the form of hydrocode snapshots. The calculation may proceed either forward in time (from bottom to top) or backwards in time (from top to bottom).
Galaxies 11 00110 g0a2
Figure A3. Simultaneous advance in two-dimensional space and forward in time of a few lines-of-sight. Top half depicts the spatial situation at t = 1 . Sixteen jet matter portions occupy this mini 4 by 4 grid. Each piece of matter is named after its position at t = 1 and retains that name as it moves along. Bottom half shows how the situation evolves as time marches on, with light rays meeting different jet segments that cross their path. A dash means a light ray meeting jet matter other than the above or nothing at all.
Figure A3. Simultaneous advance in two-dimensional space and forward in time of a few lines-of-sight. Top half depicts the spatial situation at t = 1 . Sixteen jet matter portions occupy this mini 4 by 4 grid. Each piece of matter is named after its position at t = 1 and retains that name as it moves along. Bottom half shows how the situation evolves as time marches on, with light rays meeting different jet segments that cross their path. A dash means a light ray meeting jet matter other than the above or nothing at all.
Galaxies 11 00110 g0a3
We then proceed to calculate the current length of the LOS
l los ( current ) = [ ( dlr ( nx 1 current nx 10 ) ) 2 + ( dlu ( ny 1 current ny 10 ) ) 2 + ( dlc ( nz 1 current nz 10 ) ) 2 ] 1 / 2
where LOS length is measured in cell length units and
nx 1 current = nx 10 + rc , ny 1 current = ny 10 + uc , nz 1 current = nz 10 + cc
along the x, y, and z directions, dlc, dlu, and dlr are the respective normalized hydrocode Cartesian cell lengths. Their values are usually unity or close to unity, as set in the hydrocode by the user, and rlos requires them fixed, meaning only homogeneous grids are currently supported. Furthermore, if the hydrocode grid is read by pload at a reduced resolution, rlos cell sizes are automatically adjusted accordingly.
We can finally write
l los ( current ) = [ ( ( dlr rc ) 2 ) + ( ( dlu uc ) 2 ) + ( ( dlc cc ) 2 ) ] 1 / 2
We then proceed to calculate curtime, the current hydrosimulation time of the light ray along the LOS. For forward in time LOS advance
curtime = l los ( current ) / c light + t shotmin .
while for back in time ray-tracing
curtime = + t shotmax l los ( current ) / c light .
t shotmin is the timestamp of the first loaded snapshot, t shotmax the one of the last snapshot loaded, and clight is the speed of light in cells per simulation second.
When curtime crosses a new snapshot’s time tag, the algorithm switches to drawing the LOS through the 3D volume of the new snapshot (Figure A2). We keep moving along the same LOS in 3D space, but we switch to a new time instant in the data. The above temporal shift is repeated as many times as required by the relevant criterion along the LOS until the spatial end of the LOS.

Appendix C.2.3. Aiming at the Line-of-Sight

The direction of an LOS in 3D space is defined by the two angles of azimuth (angle 1) and elevation (angle 2) (Figure 1), where the plane of angle 1 is the x y plane, parallel to xy.
LOS’s may be either parallel to each other, or focused. With a focused beam, the image is formed on a adjustable size fiducial screen located at a user-defined position between the focal point and the model system (Figure A4).
Figure A4. Geometry of focused-beam imaging in rlos, for the case of forward in time calculation. A fiducial screen is where the image is formed.
Figure A4. Geometry of focused-beam imaging in rlos, for the case of forward in time calculation. A fiducial screen is where the image is formed.
Galaxies 11 00110 g0a4
Figure A5. Geometric arrangement with regard to the viewing angles in the model for the special cases of angle2 = 0 (left) and angle1 = 0 (right). For each subcase, the arrow shows the LOS direction, which is different than the reader’s direction of view.
Figure A5. Geometric arrangement with regard to the viewing angles in the model for the special cases of angle2 = 0 (left) and angle1 = 0 (right). For each subcase, the arrow shows the LOS direction, which is different than the reader’s direction of view.
Galaxies 11 00110 g0a5
For a jet parallel to the y axis, the angle between local jet matter velocity u and LOS, losu = ( LOS , u ^ ) , is generally small when angle 1 (xy azimuth) approaches 90 degrees and vice versa (Figure A5). As is well-known [18], the angle losu affects the relativistic emission calculations. Individual jet elements may still move in directions different than the main jet axis, as part of a dynamic flow.
For a jet along the y axis, the plane of angle 2 (elevation) is largely perpendicular to the jet when azimuth (angle 1) is zero, while it is roughly parallel to the jet when azimuth (angle 1) is 90 degrees. Usually, the jet bears an approximate cylindrical symmetry, and this has an interesting effect on the sensitivity of the synthetic image to the viewing angles (Figure A5). More specifically, for a small azimuth (angle 1), if we vary elevation (angle 2), we indeed rotate the viewing point around the jet axis, thus producing similar intensities, thanks to the approximate cylindrical symmetry of the jet. Thus, for a jet moving along the y axis, the smaller azimuth (angle 1) is, the less difference varying elevation (angle 2) makes.
On the other hand, for azimuth (angle 1) nearing π /2, varying elevation (angle 2) rotates the view within a plane approximately parallel to the jet, resulting to considerable differences in the image (no symmetry involved this time). Consequently, the larger azimuth (angle 1) is, the stronger the effect, on the synthetic image, from changing elevation (angle 2).

Appendix C.3. Relativistic Effects-Doppler Boosting

Appendix C.3.1. General

The main effects of the Lorentz/Poincaré transform on the emission from a relativistic object [29], specifically applied to an astrophysical jet, are relativistic aberration, time dilation, and frequency shift [18,32,61,62].
This area refers to E/M emission, but not to particle emission. Jet spectrum is given by any suitable form inserted into the model, including the spectrum resulting from synchrotron emission and absorption coefficients. Earth frame jet spectrum S obs can be expressed [18,32] as
S obs = S jet D 3 + α
where α is the spectral index and D is the Doppler factor. Exponent (3 + α ) in the above can be broken down into different contributions from separate effects. Two units come from the aberration of light, one from the relativistic dilation of time and α from the effect of frequency shift, while for a continuous optically thin jet, a D factor is lost [32].

Aberration-Searchlight Effect

Relativistic aberration changes the perceived direction of light (there is no ray curving in special relativity) when transforming between the jet frame and the earth frame, ‘tilting rays’ emanating from the jet, generally towards its head area.
At an individual cell level, cell emission along a ray within the cell’s boost cone is then accordingly reinforced; if outside the cone, it is weakened. Depending on local velocity value and direction, neighbouring cells may have totally different boost cones.

Time Dilation

Time dilation contributes one D factor to the emission result. Again, this refers to E/M radiation.

Frequency Shift

E/M radiation emitted at a given frequency, from fast-moving jet matter, is taken to be Doppler shifted in frequency
f obs = f calc D
where f obs is the observed frequency, and f calc is the frequency used in emission calculations performed in the jet frame of reference [18].
For D ≥ 1, emission is calculated at a frequency lower than the observed, resulting in higher intensity, when the employed spectrum decreases with frequency. For radio emission, we therefore calculate emission and absorption at the f calc = f obs /D. Each computational cell, in general, has its own Doppler factor D, and therefore its own frequency shift.

Appendix C.3.2. Lorentz Factor

The Lorentz factor for a hydrocode cell is [18]
Γ Lorentz = 1 1 u 2
where
u = u x 2 + u y 2 + u z 2 1
is the value of local velocity u = ( u x , u y , u z ) , in units of the speed of light.

Appendix C.3.3. Doppler Factor Calculation

Jet radiation is either boosted or deboosted, depending on the angle ‘losu’ between the direction of the LOS and u . The higher the jet speed is, the narrower and stronger the cell boost cones are, around the direction of local velocity. On the other hand, outside cell boost cones, deboosting occurs, that is to say, the higher the velocity is, the weaker the signal becomes. D equals
D = 1 u 2 ( 1 u cos ( losu ) )
For the above, the angle between the LOS and the local velocity vector is required at every point of the computational space.
(As a note, for particles, their distribution is transformed to the Earth frame, as shown in [38]).
The cosine of angle losu is calculated in the following manner:
Let us define a fiducial unitary LOS vector ( LOS ) = ( lx 1 , lx 2 , lx 3 ) , with ( LOS ) = lx 1 2 + lx 2 2 + lx 3 2 = 1 . In the following, ϕ 1 and ϕ 2 represent azimuth and elevation angles 1 and 2, respectively.
lx 1 = cos ( ϕ 1 ) cos ( ϕ 2 ) , lx 2 = sin ( ϕ 1 ) cos ( ϕ 2 ) , lx 3 = sin ( ϕ 2 )
LOS u = ( LOS ) u cos ( LOS , u ^ ) = lx 1 u x + lx 2 u y + lx 3 u z
Therefore, we have ((LOS) = 1)
cos ( LOS , u ^ ) = lx 1 u x + lx 2 u y + lx 3 u z ( LOS ) u = lx 1 u x + lx 2 u y + lx 3 u z ( u x 2 + u y 2 + u z 2 )
For back in time ray-tracing, a minus sign is introduced to the above equation. Furthermore, a minuscule number is added to the denominator of Equation (A27), in case u = 0. The above calculation allows the assignment of a Doppler boosting factor through Equations (A23), (A24) and (A27) to each discrete emission event along a line-of-sight.

Appendix C.3.4. CoslosB Calculation

The calculation of the angle losb between the magnetic field and the LOS is performed, in rlos, in exactly the same way as for losu, above. Only this time the vector of velocity is replaced by the magnetic field one.

Alternative Frequency Shift

rlos may include different emission dependencies on frequency, where we calculate intensity at f calc and observe that at f obs .

Appendix C.4. Testing Parameters

Certain parameters, that facilitate testing rlos, are presented here.

Appendix C.4.1. The Clight Parameter

Let us consider a 4D array, comprising a succession of hydrocode snapshots. The LOS traversing those data, moves at a speed of clight cells per time unit. When we artificially adjust clight to a lower value [60,63], then the algorithm jumps to a new snapshot after spatially advancing through fewer cells. A slower LOS advances farther in time while crossing a given distance through the jet, allowing for a detailed study of the time-jumping algorithm. On the other hand, setting clight to a very high value leads to a single shot image, as we never advance to a further temporal slice.
A representative hydrocode scaling is the following (for neutrinos)
L sim = 10 10 cm , u sim = 3 × 10 10 cm s , ρ sim = 1.67 × 10 24 g cm 3 t sim = L sim u sim = 1 3 s
where t sim is the hydrocode time unit, u sim is the speed of light and L sim is the hydrocode length unit. When preparing the hydrocode run, the time span, in simulation seconds, between data snapshots, should optimally be set, to l LOS /(n*clight). l is the LOS length, in cells and n is the desired number of snapshots to cover the imaged timespan. If we employ the parameter sfactor, pload’s shrink factor, cells are enlarged and the calculated value of clight shrinks accordingly (sfactor regrids the hydrodata to a coarser grid). Overall accuracy then suffers somewhat, and regridding to a lower resolution should be used only as a preview.
Altering clight only affects the light ray speed, not the speed of matter. Consequently, overriding clight does not affect the relativistic emission calculations (like tweakspeed does, Appendix C.4.4). An altered clight is merely an artifice, introduced in post processing, in order to explore the effect of using more, or less, temporal slices in the final image.
Figure A6. A simplified flow diagram depicting the basic logical structure of rlos imaging code, for the case of rays parallel to each other. The synthetic image’s xy loops here correspond to either the yz or the xz side plane of the computational box.
Figure A6. A simplified flow diagram depicting the basic logical structure of rlos imaging code, for the case of rays parallel to each other. The synthetic image’s xy loops here correspond to either the yz or the xz side plane of the computational box.
Galaxies 11 00110 g0a6

Appendix C.4.2. The FS switch

After the hydrodata are loaded, a global operation calculates, for each cell, a jet frame frequency f calc : f calc = f obs / D , where f obs is the observing frequency and D is the local Doppler factor of the cell. The frequency shift (FS) switch selects between using the local f calc or the global f obs in the emission calculations. The FS facility allows for a direct user-defined emission dependence on frequency to be introduced ([32] pg. 199: Simplified synchrotron jet), in the form of of a function S obs = S obs ( f ) .

Appendix C.4.3. The DB Switch

The DB switch offers the option of using the Doppler boosting effect.

Appendix C.4.4. The Speed Tweak Parameter

A test is introduced, whereby matter velocity is multiplied, on a global scale, by a ‘speed tweak’ factor. This offers a quick way to observe the impact, on the synthetic image, of altering the hydrodynamic speed in post-processing, for the same simulation run. The natural value of tweakspeed is 1. At low tweak speed factors (less than 1) the effects, on the final image, of both DB and FS, are reduced, and vice versa. The maximum for tweakspeed is c/u ( max ) , above which velocites higher than c are artificially created in the grid.

Appendix C.5. rlos210

rlos version 2.10 includes a unified, functionalized, modular approach. The XZ and YZ versions were merged for both focused beam and parallel LOSs. Code was reorganised, and a series of tests are included.

Appendix C.6. rlos210 Commentary Transcript

This is the latest version of rlos. Version 2 is a major upgrade of original rlos code. This time, the programme was broken up into procedures and functions with a modular structure.
The programme allows the user to select which case to simulate, through an external parameter file. There is a unified approach, where the same modules operate on different geometries, through parameterization.
The user may select the values of the parameters of rlos version 1, and fully employ them. As mentioned above, there is no more a different version of rlos for XZ and YZ plane image formation. Now, there is one version of the code for both cases. Furthermore, for each of those cases the user may select either radiograph or camera obscura imaging technique.
The radiograph setup has all lines of sight parallel to each other, just like rlos v.1. This means the film (fiducial imaging screen) is the size of the scene (grid), like an X-ray medical image. The latter type of image shows clearly the various details of the system.
On the other hand, camera obscura, or focused beam, has a focal point where the eye of the fiducial observer is located. The imaging screen, in camera obscura, is of varied size: It may be equal, or smaller to the grid slice, at a given point along either the x or y axis depending on YZ or XZ imaging plane case. At the moment, the fiducial imaging screen must be parallel to the corresponding side of the grid, i.e., either XZ or YZ. Screen location on-axis may vary within the grid. The smaller the screen is, the smaller the image.
The focal point may reside either on the side of the grid or outside the grid but within the limits of the projection of the XZ or YZ plane. It may have a negative or zero axis position, but its two planar coordinates must be smaller than the grid size.
Direction angles are no longer necessarily constant throughout the calculation: for the focused-beam case, each LOS is drawn with a different set of azimuth (phi1) and elevation (phi2) angles. Angles are calculated using the lines that connect the focal point and the imaging screen point, which is the target point for the LOS.
The LOS then begins from the focal point if it resides on the grid side or from the LOS entry point, calculated suitably. From then on, it advances using aiming algorithms, trying to pass through the targeted screen point. It normally obtains the target or closely misses it. In general, the higher the resolution is, the better the accuracy in this respect.
For GR pseudo-Newtonian simulations, a logical next step is to introduce D(phi1), D(phi2), i.e., alter angles along an LOS from cell to cell according to the effect of the potential.
Then, innermost jet workings may be imaged, if the hydrocode can employ influence from a black hole.

Appendix C.6.1. Back in Time Integration along the LOS

In this version, calculations may be performed either ahead in time or backwards in time from a selected time instant (tpicked) backwards. For camera obscura, back in time is generally the correct way to proceed. For radiographs, ahead in time also works fine, assuming a suitable fiducial setup of the jet system vs the observer.
Tpicked is only employed when back in time switch is activated in the external parameter file. tpicked must be generally towards the end of the preselected range of dump files, or timeshots, to be loaded to RAM. At the moment, it is set equal to tmax, for convenience. A sufficient backwards time range must be provided for the LOS, in order to travel back through time without reaching the beginning of the grid. Otherwise, code cannot finish integration along the LOS. When testing, the facility of altering light speed, clight, may be used to study this effect.
Pathfinding algorithms were upgraded for this version. For each combination of XZ or YZ and radiograph or camera obscura, a certain set of such pathfinders were employed.

Appendix D. Neutrino Emission Calculations

The presentation in this section is based on the corresponding section of [23].

Appendix D.1. Proton Energy Loss

A number of energy loss mechanisms are included for the hot proton distribution [6,9,20]. The discussion of this subsection is based on a cell with the following properties: ( u x , u y , u z , b x , b y , b z , n p , ϕ 1 , ϕ 2 , α ) = (−0.3780c, 0.4480c, 0.0124c, 10 5 G, 10 6 G, 10 5 G, 2.1 × 10 11 cm 3 , 1.047 rad, 5.00 × 10 7 rad, 2.0), where u stands for velocity and b for magnetic field along directions x, y, or z. In the following, n is the jet proton density, and α is the high-energy proton distribution spectral index. As an exception, Figure A11 uses a velocity vector of (0.2, 0.8, 0.1)c. Figure A7, Figure A8, Figure A9, Figure A10, Figure A11 and Figure A12 are taken from [23].
Figure A7. Inelastic proton–proton collision standard plotted with energy. It demonstrates rather small variation (linear vertical scale) of its value over a large energy range (logarithmic horizontal scale), covering and exceeding the energy span required for the calculations that follow later in this paper.
Figure A7. Inelastic proton–proton collision standard plotted with energy. It demonstrates rather small variation (linear vertical scale) of its value over a large energy range (logarithmic horizontal scale), covering and exceeding the energy span required for the calculations that follow later in this paper.
Galaxies 11 00110 g0a7
A hot proton upper energy cutoff at E ≤ 10 6 GeV is used. The adiabatic system expansion time scale is taken to be [6]
t adb 1 = 2 3 u b ( adb ) z j
with z j = 10 11 cm being the characteristic width of the jet system at the neutrino emission size scale. For this calculation, u b ( adb ) is preset to 0.8, for a 0.8c jet.
For the p–p collision energy loss mechanism,
t pp 1 = n c σ inel pp ( E p ) K pp
where n is the bulk flow number density, K pp = 0.5 [6]. Equation (A30) is justified if we consider a small cube of matter of number density n, moving at speed (near) c, having a surface A perpendicular to its direction of motion. Then, n*c has the dimensions of cm 2 × s 1 . This is then multiplied by the inelastic p–p collision cross-section [6]
σ pp ( inel ) = ( 34.3 + 1.88 L + 0.25 L 2 ) × [ 1 ( E th E p ) 4 ] 2 × 10 27 cm 2
L = ln(E p /1000GeV), E th = 1.2 GeV. In Figure A7, σ pp ( inel ) is plotted.
The pion decay timescale is
t π = t π 0 Γ π + t esc
where Γ π is the pion Lorentz factor and
t π 0 = 2.6 × 10 8   s
Figure A8. Non-thermal proton distribution energy loss time scales, for various loss mechanisms, drawn with energy in GeV. t accel is the proton acceleration time scale at shocks. t synfvarmag stands for the synchrotron mechanism loss time scale, using a magnetic field that varies from point to point within the jet. t adb is the adiabatic loss time scale, t pp is the (hot–cold) proton–proton collision timescale. t pipion stands for the pion decay timescale t π .
Figure A8. Non-thermal proton distribution energy loss time scales, for various loss mechanisms, drawn with energy in GeV. t accel is the proton acceleration time scale at shocks. t synfvarmag stands for the synchrotron mechanism loss time scale, using a magnetic field that varies from point to point within the jet. t adb is the adiabatic loss time scale, t pp is the (hot–cold) proton–proton collision timescale. t pipion stands for the pion decay timescale t π .
Galaxies 11 00110 g0a8
In the calculations, we employ the form
t π = t π 0 ( E π m π c 2 ) + t esc
where m π is the pion mass. The light escape time t esc is a sensitive parameter of the model.
The synchrotron hot proton energy loss time scale is [6]
t sync 1 = 4 3 m e m p 3 1 8 π c m e σ T B 2 E p m p c 2
m e is the electron mass and m p the proton mass. σ T = 8 π 3 ( e 2 m e c 2 ) 2 = 6.65 × 10 25 cm 2 is the Thompson cross-section, and B is the local magnetic field. The proton’s Lorentz factor Γ p is written as E p m p c 2 , in order to facilitate energy-dependent calculations later on. In summary,
t loss 1 = t sync 1 + t adb 1 + t pp 1
The time-scales of the different energy-loss mechanisms are presented in Figure A8.
Figure A9. Density of nonthermal protons in the jet using a high-energy cutoff feature plotted with energy.
Figure A9. Density of nonthermal protons in the jet using a high-energy cutoff feature plotted with energy.
Galaxies 11 00110 g0a9
Figure A10. F π ( pp ) x , E x function, Equation (A42), corresponding to the pion spectrum emerging from a single (hot–cold) proton collision, multiplied by x= E π E p . Calculated at three different energies of the hot proton.
Figure A10. F π ( pp ) x , E x function, Equation (A42), corresponding to the pion spectrum emerging from a single (hot–cold) proton collision, multiplied by x= E π E p . Calculated at three different energies of the hot proton.
Galaxies 11 00110 g0a10
Figure A11. Pion injection function Q, weighted by pion energy, measured in non-normalized units, describing the spectrum emerging from many (non-thermal–thermal) p–p collisions. Contributions rapidly decline as particle energy increases. This figure uses a test velocity vector of (0.2, 0.8, 0.1)c.
Figure A11. Pion injection function Q, weighted by pion energy, measured in non-normalized units, describing the spectrum emerging from many (non-thermal–thermal) p–p collisions. Contributions rapidly decline as particle energy increases. This figure uses a test velocity vector of (0.2, 0.8, 0.1)c.
Galaxies 11 00110 g0a11

Appendix D.2. Particle Cascades in the Jets

Hot–cold proton interaction results to a distribution of high-energy pions, which then decay, allowing for the creation of energetic neutrinos. We have [43,47,64,65]
p p p p π 0 + π 0 ,
for neutral pions π 0 , and
p p p n π + , p p p n π + π + + π + ,
for π ± .
Neutral pions π 0 decay to gamma rays. On the other hand, π ± mainly decay to an antimuon and a muonic neutrino, or to a muon and an antineutrino (prompt neutrinos) [47,65]. We therefore obtain neutrinos that escape the system, after the cascade.
π + μ + + ν μ , π μ + ν ˜ μ .
As an approximation, we do not consider neutrino production through secondary channels or delayed neutrinos.
For each successive particle population in the above cascades, the transport equation for nonstochastic phenomena and for time-independent transport (transport time much less than the time step of the dynamic simulation), takes a simplified form:
N E + N t loss = Q ( E , r )
where r is the spatial position vector, N is the particle density of the daughter population, and Q is its injection function.
A power-law distribution is assumed for protons (Figure A9 shows power-law fast proton density), replacing the solution of the first transport equation in the cascade. Afterwards, along the cascade, we calculate the properties of resulting particle distributions [20].
Figure A12. Pion energy distribution plotted in non-normalized units versus energy. In the software, the above distribution is represented by function U.
Figure A12. Pion energy distribution plotted in non-normalized units versus energy. In the software, the above distribution is represented by function U.
Galaxies 11 00110 g0a12

Appendix D.3. Lorentz Transform of High E Proton Distribution

For the calculation of the fast proton distribution, the relevant directional equation (direction is defined by the angle θ between velocity and line of sight) can be found in [38,66]. The selected variant originates from [38], minus a geometry factor that is here absorbed into the normalization factor [23]
n ( E , θ ) = Γ α 1 E α ( 1 β cos ( θ ) 1 m 2 c 4 E 2 ) α 1 [ sin 2 ( θ ) + Γ 2 ( cos ( θ ) β 1 m 2 c 4 E 2 ) 2 ] 1 2
where Γ is the Lorentz factor of the particles in a tiny volume.

Appendix D.4. Pion Injection Function and Pion Energy Distribution

For each fast–slow proton interaction, a spectrum of possible pion energies exists, given by function F π [6,9,20].
F π ( pp ) x , E x = 4 α B π x α 1 1 x α 1 + r x α ( 1 x α ) 4 1 1 x α + r ( 1 2 x α ) 1 + r x α ( 1 x α ) 1 m π c 2 x E p 1 2
where x = E / E p , B π = α + 0.25, α = 3.67 +0.83L + 0.075L 2 , r = 2.6/ α , α = 0.98/ α [6,20].
Figure A10 x F π is plotted with the fraction x for different fast proton energies.
Pion injection function Q π ( pp ) comprises pion contributions at each pion energy to that pion energy from spectrum F of all potential p–p interactions.
Q π ( pp ) ( E , z ) = n ( z ) c E E p ( max ) 1 d x x E x , z F π ( pp ) x , E x σ pp ( inel ) E x ,
x is the fraction of the pion energy to proton energy, and n ( z ) is the jet flow proton density.
Figure A11 plots Q π ( pp ) versus pion energy E π .
In order to obtain the pion distribution, the following transport equation is employed:
N π E + N π t loss = Q π ( pp ) ( E , z )
where N π ( E , z ) denotes the pion energy distribution. Then
N π ( E ) = 1 | b π ( E ) | E E ( max ) d E Q π ( pp ) ( E ) exp [ τ π ( E , E ) ] ,
and
τ π ( E , E ) = E E d E t π 1 ( E ) | b π ( E ) | .
b π ( E ) = E ( t sync 1 + t adb 1 + t π p 1 + t π γ 1 ) is the energy loss rate of the pion. As an approximation, the last term in the latter expression is omitted. In Figure A12 plots nemiss [49] software function U (U analytical ), representing N π , with pion energy.
The above calculations are performed for each computational cell. A cell is macroscopically large inasmuch as only the deterministic portion of the transport equation is employed, in turn rendering it deterministic. Again, we 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 taken to be smaller enough than the hydrocode’s time step, so that a cell’s radiative interactions belong to a single time step each time.

Appendix D.5. Neutrino Emissivity

As mentioned in the main text, the emissivity of prompt neutrinos [5,6,20,21], is
Q π ν ( E ) = E E max d E π t π 1 ( E π ) N π ( E π ) Θ ( 1 r π x ) E π ( 1 r π ) ,
E is neutrino energy, x = E / E π , r π =(m μ /m π ) 2 and t π is the pion decay timescale. Θ ( χ ) is the theta function [6,43]. Neutrino emissivity is calculated for each individual cell using the angle of the local velocity to the LOS crossing that cell.

Appendix E. Flow Diagrams

In this section, some code flow diagrams are provided, as supplemental material.
Figure A13. Flow diagram of nemiss depicting the basic structure of the programme.
Figure A13. Flow diagram of nemiss depicting the basic structure of the programme.
Galaxies 11 00110 g0a13
Figure A14. Flow diagram of the main loop of nemiss where the calculation takes place. This was upgraded to 5D, with the addition of a loop for particle energy.
Figure A14. Flow diagram of the main loop of nemiss where the calculation takes place. This was upgraded to 5D, with the addition of a loop for particle energy.
Galaxies 11 00110 g0a14
Figure A15. Flow diagram of rlos.1 depicting procedures and functions called during programme execution.
Figure A15. Flow diagram of rlos.1 depicting procedures and functions called during programme execution.
Galaxies 11 00110 g0a15

References

  1. Mirabel, I.F.; Rodríguez, L.F. Sources of relativistic jets in the Galaxy. Ann. Rev. Astron. Astrophys. 1999, 37, 409. [Google Scholar] [CrossRef]
  2. Romero, G.E.; Torres, D.F.; Kaufman Bernadó, M.M.; Mirabel, I.F. Hadronic gamma-ray emission from windy microquasars. Astron. Astrophys. 2003, 410, L1–L4. [Google Scholar] [CrossRef]
  3. Bednarek, W. TeV Neutrinos from Microquasars in Compact Massive Binaries. Astrophys. J. 2005, 631, 466–470. [Google Scholar] [CrossRef]
  4. Bosch-Ramon, V. Theoretical overview on high-energy emission in microquasars. Astrophys. Space Sci. 2007, 309, 321–331. [Google Scholar] [CrossRef]
  5. Reynoso, M.M.; Romero, G.E.; Christiansen, H.R. Production of gamma rays and neutrinos in the dark jets of the microquasar SS433. Mon. Not. R. Astron. Soc. 2008, 387, 1745–1754. [Google Scholar] [CrossRef]
  6. Reynoso, M.M.; Romero, G.E. Magnetic field effects on neutrino production in microquasars. Astron. Astrophys. 2009, 493, 1–11. [Google Scholar] [CrossRef]
  7. Christiansen, H.R. High energy emission from galactic jets. arXiv 2013, arXiv:1306.1792. [Google Scholar]
  8. Zhang, J.F.; Feng, Y.G.; Lei, M.C.; Tang, Y.Y.; Tian, Y.P. High-energy neutrino emission from low-mass microquasars. Mon. Not. R. Astron. Soc. 2010, 407, 2468–2474. [Google Scholar] [CrossRef]
  9. Reynoso, M.M.; Carulli, A.M. On the possibilities of high-energy neutrino production in the jets of microquasar SS433 in light of new observational data. Astropart. Phys. 2019, 109, 25–32. [Google Scholar] [CrossRef]
  10. Carulli, A.M.; Reynoso, M.M.; Romero, G.E. Neutrino production in Population III microquasars. Astropart. Phys. 2021, 128, 102557. [Google Scholar] [CrossRef]
  11. Escobar, G.J.; Pellizza, L.J.; Romero, G.E. Highly collimated microquasar jets as efficient cosmic-ray sources. Astron. Astrophys. 2022, 665, A145. [Google Scholar] [CrossRef]
  12. Gutiérrez, E.M.; Combi, L.; Romero, G.E.; Campanelli, M. Flares from dual jets in supermassive black hole binaries. arXiv 2023, arXiv:2301.04280. [Google Scholar]
  13. Egron, E.; Rodriguez, J.; Trushkin, S.A.; Pellizzoni, A.; Pilia, M.; Melis, A.; Righini, S.; Bouchet, T.; Cangemi, F.; Coleiro, A.; et al. A bright radio flare observed in GRS 1915+105 with the Medicina and RATAN radio telescopes with no X-ray counterpart in the INTEGRAL energy range. Astron. Telegr. 2023, 16008, 1. [Google Scholar]
  14. Koessl, D.; Mueller, E.; Hillebrandt, W. Numerical simulations of axially symmetric magnetized jets. I—The influence of equipartition magnetic fields. II—Apparent field structure and theoretical radio maps. III—Collimation of underexpanded jets by magnetic fields. Astron. Astrophys. 1990, 229, 378–415. [Google Scholar]
  15. Rieger, F.M.; Duffy, P. A Microscopic Analysis of Shear Acceleration. Astrophys. J. 2006, 652, 1044–1049. [Google Scholar] [CrossRef]
  16. Rieger, F.M. An Introduction to Particle Acceleration in Shearing Flows. Galaxies 2019, 7, 78. [Google Scholar] [CrossRef]
  17. Kuldeep, S. Study of relativistic magnetized outflows with relativistic equation of state. Mon. Not. R. Astron. Soc. 2019, 488, 5713–5727. [Google Scholar] [CrossRef]
  18. Hughes, P.A. Beams and Jets in Astrophysics; Cambridge University Press: Cambridge, UK, 1991. [Google Scholar]
  19. Mignone, A.; Bodo, G.; Massaglia, S.; Matsakos, T.; Tesileanu, O.; Zanni, C.; Ferrari, A. PLUTO: A numerical code for computational astrophysics. Astrophys. J. Supp. 2007, 170, 228. [Google Scholar] [CrossRef]
  20. Kelner, S.R.; Aharonian, F.A.; Bugayov, V.V. Energy spectra of gamma rays, electrons, and neutrinos produced at proton-proton interactions in the very high energy regime. Phys. Rev. D 2006, 74, 034018. [Google Scholar] [CrossRef]
  21. Lipari, P.; Lusignoli, M.; Meloni, D. Flavor composition and energy spectrum of astrophysical neutrinos. Phys. Rev. D 2007, 75. [Google Scholar] [CrossRef]
  22. Pacholczyk, A.G. Radio Astrophysics. Nonthermal Processes in Galactic and Extragalactic Sources. Phys. Today 1971, 24, 57–58. [Google Scholar] [CrossRef]
  23. Smponias, T. Synthetic Neutrino Imaging of a Microquasar. Galaxies 2021, 9, 80. [Google Scholar] [CrossRef]
  24. Lampa, A. Wie erscheint nach der Relativitätstheorie ein bewegter Stab einem ruhenden Beobachter? Z. Phys. 1924, 27, 138. [Google Scholar] [CrossRef]
  25. Terrell, J. Invisibility of the Lorentz Contraction. Phys. Rev. 1959, 116, 1041. [Google Scholar] [CrossRef]
  26. Penrose, R. The apparent shape of a relativistically moving sphere. Proc. Camb. Phil. Soc. 1959, 55, 137. [Google Scholar] [CrossRef]
  27. Weisskopf, V.F. The visual appearance of rapidly moving objects. Phys. Today 1960, 13, 24. [Google Scholar]
  28. Scott, G.D.; Viner, M.R. The Geometrical Appearance of Large Objects Moving at Relativistic Speeds. Am. J. Phys. 1965, 33, 534. [Google Scholar] [CrossRef]
  29. Weiskopf, D. Visualization of Four-Dimensional Spacetimes. Ph.D. Thesis, der Eberhard-Karls-Universitat zu Tubingen, der Fakultat fur Physik, Tübingen, Germany, 2001. [Google Scholar]
  30. Deissler, R.J. The appearance, apparent speed, and removal of optical effects for relativistically moving objects. Am. J. Phys. 2005, 73, 663. [Google Scholar] [CrossRef]
  31. Weiskopf, D. Scientific Visualization: Advanced Concepts; Hagen, H., Ed.; Dagstuhl Publishing, Leibniz Center for Informatics: Wadern, Germany, 2010; pp. 289–302. [Google Scholar]
  32. Cawthorne, T.V. Interpretation of Parsec scale jets. In Beams and Jets in Astrophysics; Cambridge University Press: Cambridge, UK, 1991. [Google Scholar]
  33. Weiskopf, D.; Kraus, U.; Ruder, H. Searchlight and Doppler Effects in the Visualization of Special Relativity: A Corrected Derivation of the Transformation of Radiance. ACM Trans. Graph. 1999, 18, 3. [Google Scholar] [CrossRef]
  34. Klauber, R.D. Is detection of Fitzgerald-Lorentz contraction possible? arXiv 2008, arXiv:0808.1117v1. [Google Scholar]
  35. Rybicki, G.B.; Lightman, A.P. Radiative Processes in Astrophysics; Wiley-Interscience: New York, NY, USA, 1979. [Google Scholar]
  36. Jarabo, A.; Masia, B.; Velten, A.; Barsi, C.; Raskar, R.; Gutierrez, D. Relativistic Effects for Time-Resolved Light Transport. Comput. Graph. Forum 2015, 34, 1–12. [Google Scholar] [CrossRef]
  37. Hickey, F.R. Two-dimensional appearance of a relativistic cube. Am. J. Phys. 1979, 47, 711. [Google Scholar] [CrossRef]
  38. Torres, D.F.; Reimer, A. Hadronic beam models for quasars and microquasars. Astron. Astrophys. 2011, 528, L2. [Google Scholar] [CrossRef]
  39. Lobato, R.V.; Coelho, J.G.; Malheiro, M. Ultra-high energy cosmic rays from white dwarf pulsars and the Hillas criterion. J. Phys. Conf. Ser. 2017, 861, 012005. [Google Scholar] [CrossRef]
  40. Begelman, M.C.; Blandford, R.D.; Rees, M.J. Massive black hole binaries in active galactic nuclei. Nature 1980, 287, 307–309. [Google Scholar] [CrossRef]
  41. Romero, G.E.; Müller, A.L.; Roth, M. Particle acceleration in the superwinds of starburst galaxies. Astron. Astrophys. 2018, 616, A57. [Google Scholar] [CrossRef]
  42. Derishev, E.V.; Aharonian, F.A.; Kocharovsky, V.V. High-energy emission from off-axis relativistic jets. In Proceedings of the High Energy Gamma-Ray Astronomy, Heidelberg, Germany, 26–30 July 2004; Aharonian, F.A., Völk, H.J., Horns, D., Eds.; American Institute of Physics Conference Series. AIP: New York, NY, USA, 2005; Volume 745, pp. 510–515. [Google Scholar] [CrossRef]
  43. Smponias, T.; Kosmas, O. High Energy Neutrino Emission from Astrophysical Jets in the Galaxy. Adv. High Energy Phys. 2015, 2015, 921757. [Google Scholar] [CrossRef]
  44. Smponias, T. RLOS: Time-Resolved Imaging of Model Astrophysical Jets. 2018. Available online: https://ascl.net/1811.009 (accessed on 23 October 2023).
  45. Smponias, T.; Kosmas, T.S. Modelling the equatorial emission in a microquasar. Mon. Not. R. Astron. Soc. 2011, 412, 1320. [Google Scholar] [CrossRef]
  46. Smponias, T.; Kosmas, T.S. Dynamical and radiative simulations of γ-ray jets in microquasars. Mon. Not. R. Astron. Soc. 2014, 438, 1014. [Google Scholar] [CrossRef]
  47. Kosmas, O.; Smponias, T. Simulations of Gamma-Ray Emission from Magnetized Microquasar Jets. arXiv 2018, arXiv:1808.00303. [Google Scholar] [CrossRef]
  48. Hjellming, R.M.; Johnston, K.J. Radio Emission from Conical Jets Associated with X-Ray Binaries. ApJ 1988, 328, 600. [Google Scholar] [CrossRef]
  49. Smponias, T. nemiss: Neutrino Imaging of Model Astrophysical Jets. 2019. Available online: github.com/teoxxx/nemiss_pbl (accessed on 23 October 2023).
  50. Fabrika, S. The jets and supercritical accretion disk in SS433. Astrophys. Space Phys. Rev. 2004, 12, 1–152. [Google Scholar] [CrossRef]
  51. Kološ, M.; Tursunov, A.; Stuchlík, Z. Possible signature of the magnetic fields related to quasi-periodic oscillations observed in microquasars. Eur. Phys. J. C 2017, 77, 860. [Google Scholar] [CrossRef]
  52. Paredes, J.M.; Marti, J. Microquasars in the galaxy. Contrib. Sci. 2003, 2, 303–314. [Google Scholar]
  53. Ambrogi, L.; Celli, S.; Aharonian, F. On the potential of Cherenkov Telescope Arrays and KM3 Neutrino Telescopes for the detection of extended sources. Astropart. Phys. 2018, 100, 69–79. [Google Scholar] [CrossRef]
  54. Hsiung, P.; Dunn, R. Visualizing relativistic effects in space time. In Proceedings of the Supercomputing 1989 Conference SC’89: International Conference for High Performance Computing, Networking, Storage and Analysis, Reno, NV, USA, 12–17 November 1989; Volume 1, p. 597. [Google Scholar]
  55. Hsiung, P.; Thibadeau, R.; Wu, M. T-Buffer: Fast Visualization of Relativistic Effects in Spacetime. Comput. Graph. 1990, 24, 83. [Google Scholar] [CrossRef]
  56. Kraus, U. First-person visualizations of the special and general theory of relativity. Am. J. Phys. 2000, 68, 56. [Google Scholar] [CrossRef]
  57. Muller, T.; Boblest, S. Visual appearance of wireframe objects in special relativity. Eur. J. Phys. 2014, 35, 06502. [Google Scholar] [CrossRef]
  58. Muller, T. GeoViS—Relativistic ray tracing in four-dimensional spacetimes. Comput. Phys. Commun. 2014, 185, 2301. [Google Scholar] [CrossRef]
  59. Velten, A.; Wu, D.; Jarabo, A.; Masia, B.; Barsi, C.; Joshi, C.; Lawson, E.; Bawendi, M.; Gutierrez, D.; Raskar, R. Femto-photography: Capturing and visualizing the propagation of light. ACM Trans. Graph. 2013, 32, 4. [Google Scholar] [CrossRef]
  60. Sherin, Z.W.; Cheu, R.; Tan, P.; Kortemeyer, G. Visualizing Relativity: The OpenRelativity Project. Am. J. Phys. 2016, 84, 369. [Google Scholar] [CrossRef]
  61. Sparks, W.B.; Fraix-Burnet, D.; Macchetto, F.; Owen, F.N. A counterjet in the elliptical galaxy M87. Nature 1992, 355, 804. [Google Scholar] [CrossRef]
  62. Laing, R.; Bridle, A.H. Relativistic models and the jet velocity field in the radio galaxy 3C 31. Mon. Not. R. Astron. Soc. 2002, 336, 328. [Google Scholar] [CrossRef]
  63. Kortemeyer, G.; Tan, P.; Schirra, S. A Slower Speed of Light: Developing intuition about special relativity with games FDG 2013. In Proceedings of the FDG’13—International Conference on the Foundations of Digital Games, Chania, Crete, Greece, 14–17 May 2013; ACM: New York, NY, USA, 2013; Volume 1, p. 400. [Google Scholar]
  64. Smponias, T.; Kosmas, O. Neutrino Physics in the Frontiers of Intensities and Very High Sensitivities 2016. Adv. High Energy Phys. 2017, 2017, 4962741. [Google Scholar]
  65. Campion, S.; Fuksman, J.D.M.; Hernandez, J.A.R. Neutrino production from proton-proton interactions in binary-driven hypernovae. arXiv 2020, arXiv:1910.10439. [Google Scholar]
  66. Purmohammad, D.; Samimi, J. On the hadronic beam model of TeV gamma-ray flares from blazars. Astron. Astrophys. 2001, 371, 61–67. [Google Scholar] [CrossRef]
Figure 2. The radio-scale light jet model, shown at snapshot 50 (t = 100 ks). The plasmoids traverse the stellar wind, moving in opposite directions. The approaching jet is the one moving towards the beginning of the axes. The ambient wind is also proportionally lighter than in the heavy jet model, therefore resulting in similar dynamics to the heavier jet model. The plasmoids here also traverse a dense inner stellar wind, giving rise to an equatorial concentration of matter, detectable in the synthetic radio images. Overall, this lighter jet is more realistic in terms of densities and overall kinetic luminosity. In this figure, length units are taken from the MHD simulation, where two such units equal a computational cell in length. One MHD length unit here equals 10 13 cm (Table 1).
Figure 2. The radio-scale light jet model, shown at snapshot 50 (t = 100 ks). The plasmoids traverse the stellar wind, moving in opposite directions. The approaching jet is the one moving towards the beginning of the axes. The ambient wind is also proportionally lighter than in the heavy jet model, therefore resulting in similar dynamics to the heavier jet model. The plasmoids here also traverse a dense inner stellar wind, giving rise to an equatorial concentration of matter, detectable in the synthetic radio images. Overall, this lighter jet is more realistic in terms of densities and overall kinetic luminosity. In this figure, length units are taken from the MHD simulation, where two such units equal a computational cell in length. One MHD length unit here equals 10 13 cm (Table 1).
Galaxies 11 00110 g002
Figure 3. The radio-scale heavy jet model, shown at snapshot 50 (t = 100 ks). Blobs traverse the stellar wind, moving in opposite directions and forming a pair of intermittent jets. The approaching jet is the one moving towards the beginning of the axes. A cocoon formation appears, inflated by fast moving bolides, traversing the ambient wind of the companion star, which is denser than the jet near its base. Inner wind material is pushed sideways by the jet, facilitating the creation of an equatorial zone, also detected in the synthetic radio images of the model system. In this figure, MHD simulation length units are employed, where two such units equal a computational cell in length. An MHD length unit here equals 10 13 cm (Table 1).
Figure 3. The radio-scale heavy jet model, shown at snapshot 50 (t = 100 ks). Blobs traverse the stellar wind, moving in opposite directions and forming a pair of intermittent jets. The approaching jet is the one moving towards the beginning of the axes. A cocoon formation appears, inflated by fast moving bolides, traversing the ambient wind of the companion star, which is denser than the jet near its base. Inner wind material is pushed sideways by the jet, facilitating the creation of an equatorial zone, also detected in the synthetic radio images of the model system. In this figure, MHD simulation length units are employed, where two such units equal a computational cell in length. An MHD length unit here equals 10 13 cm (Table 1).
Galaxies 11 00110 g003
Figure 4. Unnormalized synthetic radio synchrotron images of the model system at 8 GHz. Top: the heavier jet model at snapshot 60 (t = 120 ks). Bottom: the lighter jet at snapshot 105 (t = 210 ks). In both cases, we can observe blobs moving in diametrically opposite directions. A finite c is taken into consideration in these plots, resulting in the approaching blob apparently moving much faster than the receding one. The approaching plasmoid is also much brighter, even though the blobs in each pair are essentially the same in the hydrocode. Furthermore, the images shown here correspond to earlier blob locations in the hydrocode run due to the delay in the ray arrival to the fiducial observer. The total intensity in the heavy model is roughly two orders of magnitude larger than the lighter one, in rough proportion to the ratio of jet nozzle densities in the two cases. Total intensity for each image appears at the top left hand corner. A stronger equatorial emission appears in the heavy jet case, attributed to the higher densities of this run, as well as to the earlier timetag of the heavy model synthetic image.
Figure 4. Unnormalized synthetic radio synchrotron images of the model system at 8 GHz. Top: the heavier jet model at snapshot 60 (t = 120 ks). Bottom: the lighter jet at snapshot 105 (t = 210 ks). In both cases, we can observe blobs moving in diametrically opposite directions. A finite c is taken into consideration in these plots, resulting in the approaching blob apparently moving much faster than the receding one. The approaching plasmoid is also much brighter, even though the blobs in each pair are essentially the same in the hydrocode. Furthermore, the images shown here correspond to earlier blob locations in the hydrocode run due to the delay in the ray arrival to the fiducial observer. The total intensity in the heavy model is roughly two orders of magnitude larger than the lighter one, in rough proportion to the ratio of jet nozzle densities in the two cases. Total intensity for each image appears at the top left hand corner. A stronger equatorial emission appears in the heavy jet case, attributed to the higher densities of this run, as well as to the earlier timetag of the heavy model synthetic image.
Galaxies 11 00110 g004
Figure 5. Apparent superluminal motion study setup; the jet is at 55 degrees to the LOS with a speed of 0.8c. In the synthetic images (top left-hand corner and bottom left-hand corner), approaching and receding blobs are shown for a ray speed of ≃200c (top left) at shot number 30, and for normal c (bottom left) at shot number 60 (the synthetic image view is rotated around the z axis by 180 degrees for visual clarity). The difference in timing between the two synthetic images represents an attempt to match the time delay of the normal c image to the single-shot image at 200c. Distances in the synthetic images are in computational cells, where one cell equals two hydrocode length units. Top left: no real difference exists between approaching and receding plasmoids, resembling the corresponding hydrocode density plot to the right. Bottom left: the approaching blob on the fiducial screen (sky plane) appears to be much faster than the receding blob; a clear demonstration of the apparent acceleration effect. On their right-hand side, a hydrocode data rendering of the same model run is shown; the scale is in hydrocode length units. A hydrocode density plot is also shown, demonstrating the inherently symmetric motion of approaching and receding hydrocode blobs on the fiducial sky plane (without using an imaging code).
Figure 5. Apparent superluminal motion study setup; the jet is at 55 degrees to the LOS with a speed of 0.8c. In the synthetic images (top left-hand corner and bottom left-hand corner), approaching and receding blobs are shown for a ray speed of ≃200c (top left) at shot number 30, and for normal c (bottom left) at shot number 60 (the synthetic image view is rotated around the z axis by 180 degrees for visual clarity). The difference in timing between the two synthetic images represents an attempt to match the time delay of the normal c image to the single-shot image at 200c. Distances in the synthetic images are in computational cells, where one cell equals two hydrocode length units. Top left: no real difference exists between approaching and receding plasmoids, resembling the corresponding hydrocode density plot to the right. Bottom left: the approaching blob on the fiducial screen (sky plane) appears to be much faster than the receding blob; a clear demonstration of the apparent acceleration effect. On their right-hand side, a hydrocode data rendering of the same model run is shown; the scale is in hydrocode length units. A hydrocode density plot is also shown, demonstrating the inherently symmetric motion of approaching and receding hydrocode blobs on the fiducial sky plane (without using an imaging code).
Galaxies 11 00110 g005
Figure 6. Radio synchrotron model time curves, presenting the evolution of the total unnormalized intensity over the injection of a new plasmoid pair. Blob injection occurs during the first 30 out of every 100 time units. In this figure, the injection of the second pair of blobs is presented. The delay in their detection, of around 25 time units, is attributed to the finite ray travel time from the model system to the fiducial imaging plane. In both cases, a faster rise is followed by a more gradual decline. Nozzle densities: heavy jet, ρ jet = 10 12 cm 3 ; light jet, ρ jet = 10 10 cm 3 .
Figure 6. Radio synchrotron model time curves, presenting the evolution of the total unnormalized intensity over the injection of a new plasmoid pair. Blob injection occurs during the first 30 out of every 100 time units. In this figure, the injection of the second pair of blobs is presented. The delay in their detection, of around 25 time units, is attributed to the finite ray travel time from the model system to the fiducial imaging plane. In both cases, a faster rise is followed by a more gradual decline. Nozzle densities: heavy jet, ρ jet = 10 12 cm 3 ; light jet, ρ jet = 10 10 cm 3 .
Galaxies 11 00110 g006
Figure 7. Synthetic model radio images of an approaching rectangular thin plasmoid. Top: the center of the rectangle has a higher intensity that its edges, attributed to a faster arrival of the central rays compared to rays from the edges of the rectangle. The intensity peak is not fully symmetric within the inner rectangle as the imaging focal point is not exactly on the plasmoid’s trajectory. On the other hand, by setting a speed of light c to much higher than normal (bottom), this effect largely vanishes. We can also see shadows of the central region and of the receding twin jet blob. The vertical scale for intensity is in linear, arbitrary units.
Figure 7. Synthetic model radio images of an approaching rectangular thin plasmoid. Top: the center of the rectangle has a higher intensity that its edges, attributed to a faster arrival of the central rays compared to rays from the edges of the rectangle. The intensity peak is not fully symmetric within the inner rectangle as the imaging focal point is not exactly on the plasmoid’s trajectory. On the other hand, by setting a speed of light c to much higher than normal (bottom), this effect largely vanishes. We can also see shadows of the central region and of the receding twin jet blob. The vertical scale for intensity is in linear, arbitrary units.
Galaxies 11 00110 g007
Figure 8. A continuous jet is shown on the neutrino emission scale, traversing the heavier accretion disk wind construct and also the companion star’s stellar wind. Narrow collimated jets, lighter than the inner winds, propagate first through the accretion disk wind and then the stellar wind, thus inflating a “composite” cocoon structure. In this model, the accretion disk construct appears to divide the inner jet region into two separate dynamical environments, one for each jet.
Figure 8. A continuous jet is shown on the neutrino emission scale, traversing the heavier accretion disk wind construct and also the companion star’s stellar wind. Narrow collimated jets, lighter than the inner winds, propagate first through the accretion disk wind and then the stellar wind, thus inflating a “composite” cocoon structure. In this model, the accretion disk construct appears to divide the inner jet region into two separate dynamical environments, one for each jet.
Galaxies 11 00110 g008
Figure 9. Neutrino unnormalized intensity plots at snapshots 45 (t = 90 s, top) and 90 (t = 180 s, bottom). The sum of all synthetic neutrino image pixels at each energy is plotted at each time instant. At the earlier time instant, the intensities are higher across the spectrum. This is attributed to intensity decreasing with time. The finite nature of the speed of light was taken into account when calculating these plots. Furthermore, at both time instants the intensity decreases with energy over an energy spectrum covering six orders of magnitude. In this model run, jet dynamics evolve smoothly, so the spectra remain very similar in shape across the two instants here.
Figure 9. Neutrino unnormalized intensity plots at snapshots 45 (t = 90 s, top) and 90 (t = 180 s, bottom). The sum of all synthetic neutrino image pixels at each energy is plotted at each time instant. At the earlier time instant, the intensities are higher across the spectrum. This is attributed to intensity decreasing with time. The finite nature of the speed of light was taken into account when calculating these plots. Furthermore, at both time instants the intensity decreases with energy over an energy spectrum covering six orders of magnitude. In this model run, jet dynamics evolve smoothly, so the spectra remain very similar in shape across the two instants here.
Galaxies 11 00110 g009
Figure 10. Unnormalized neutrino intensity images of the same hydrocode run at snapshots 45 (top; t = 90 s) and 90 (bottom; t = 180 s) at a specific energy (400 GeV). In order to produce the image in rlos, the focused beam imaging method, with back in time integration along the LOS, was employed. In all other energy slots of the spectral energy distribution, the image shape is quite similar, though not same, but the intensity scale does fall with energy. Emission calculations were double-filtered by setting a maximum angle between the local LOS and local u and a minimum local velocity.
Figure 10. Unnormalized neutrino intensity images of the same hydrocode run at snapshots 45 (top; t = 90 s) and 90 (bottom; t = 180 s) at a specific energy (400 GeV). In order to produce the image in rlos, the focused beam imaging method, with back in time integration along the LOS, was employed. In all other energy slots of the spectral energy distribution, the image shape is quite similar, though not same, but the intensity scale does fall with energy. Emission calculations were double-filtered by setting a maximum angle between the local LOS and local u and a minimum local velocity.
Galaxies 11 00110 g010
Figure 11. The absolute value of the toroidal magnetic field component of the continuous jet is shown on the neutrino emission scale. The scale is logarithmic and the field appears to thread the jet, contributing to its confinement along its path. The field also plays an important role in the neutrino emission calculations.
Figure 11. The absolute value of the toroidal magnetic field component of the continuous jet is shown on the neutrino emission scale. The scale is logarithmic and the field appears to thread the jet, contributing to its confinement along its path. The field also plays an important role in the neutrino emission calculations.
Galaxies 11 00110 g011
Table 1. Three different imaging runs based on three separate underlying hydrocode runs.
Table 1. Three different imaging runs based on three separate underlying hydrocode runs.
ModelRadio HeavyRadio Light ν -ScaleComments
l cell ( × 10 10 cm)2.0 × 10 3 2.0  × 10 3 2.0PLUTO cell
ρ jet (cm 3 ) 1.0 × 10 12 1.0 × 10 10 1.0 × 10 11 Jet matter density
ρ w (cm 3 ) 1.0 × 10 13 1.0 × 10 11 1.0 × 10 13 Max wind density
time unit (s) 1.0 × 10 3 1.0 × 10 3 1.0 Model time scale (model s)
ρ dw (cm 3 )-- 2.0 × 10 13 Max disk wind density
t run max (s)242  × 10 3 242  × 10 3 204Model run time
MethodP. L. P. L. P. L. Piecewise linear
IntegratorM. H. M. H. M. H. MUSCL-Hancock
EOSIdealIdealIdealEquation of state
physicsRelMHDRelMHDRelMHDPLUTO setup
B field (G)1010 1.0 × 10 4 Initial toroidal magnetic field
BinSep (cm)subcellsubcell 4.0 × 10 12 Binary separation
M BH / M -5–203–10VE compact star mass
M / M -10–3010–30Companion mass
β = v 0 / c 0.80.80.8Initial jet speed
L k p 10 44 10 42 2 × 10 38 Jet kinetic luminosity
Jet typeint.int.cont.intermittent or continuous
L k ( av ) p ≃10 43 ≃10 41 2 × 10 38 Average Jet kinetic luminosity
Grid resolution60 × 100 × 5060 × 100 × 5060 × 100 × 50PLUTO grid size (xyz)
Imaging methodFBFBFBFocused beam
Time delayonononNormal ray speed
Imaging planeYZYZYZFiducial screen parallel to YZ
Emissionradio-syncradio-syncneutrinosSynthetic emission type
Code usedPLUTO-rlosPLUTO-rlosPLUTO-nemiss-rlospipeline portion employed
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Smponias, T. Simulated Radio and Neutrino Imaging of a Microquasar. Galaxies 2023, 11, 110. https://doi.org/10.3390/galaxies11060110

AMA Style

Smponias T. Simulated Radio and Neutrino Imaging of a Microquasar. Galaxies. 2023; 11(6):110. https://doi.org/10.3390/galaxies11060110

Chicago/Turabian Style

Smponias, Theodoros. 2023. "Simulated Radio and Neutrino Imaging of a Microquasar" Galaxies 11, no. 6: 110. https://doi.org/10.3390/galaxies11060110

APA Style

Smponias, T. (2023). Simulated Radio and Neutrino Imaging of a Microquasar. Galaxies, 11(6), 110. https://doi.org/10.3390/galaxies11060110

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop