Simulated Radio and Neutrino Imaging of a Microquasar

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.


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.

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.

Radio Synchrotron Emission and Self-Absorption
Adopting high-energy electron and proton acceleration at shock fronts [22], a powerlaw 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, where K pl is a power law constant for high-energy jet electrons.K pl ρ 4d = N 0 ; therefore, the reference electron density N 0 is taken as proportional to the local bulk flow proton density ρ 4d based on charge neutrality.ρ 4d 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 = ZqBR 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]: where f pp(4d) is a proton profile function, with a threshold velocity u th of 0. (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 where Γ is the Gamma function.
The above coefficients are computed at each point (computational cell) of the model jet system.The angle ( − − → LOS, − → B 4d ) 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.

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]: 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 = KN 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].
In this work, the high-energy proton distribution is considered isotropic in the jet frame.The hypothetical anisotropy of the hot proton distribution can be reflected in theneutrino 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.

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 a pion distribution and then a neutrino population, which then leave the system [5,6,20,21].The resulting neutrino population is 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.

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 Figures A1 and A4 in the Appendix C).Tracing their way along the given direction, they reach a length of (x 2 max + y 2 max + z 2 max ), 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.

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.

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.

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.

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.

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].

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 ⃗ rz = 8πρ ⃗ rz [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 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.

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.

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.

Radio-Scale Model Results
In the radio-scale model, ejected plasmoids form an intermittent twin jet (Figures 2 and 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 (Figures 2 and 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 (Figures 4 and 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 (Figures 4 and 5).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).
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).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).
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.

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 be-forehand 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 halfangle 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].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.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 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 blobson the fiducial sky plane (without using an imaging code).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.
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.
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.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 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.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].
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 neutrinoemitting 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.

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.
The intensity of the jet is then expressed as I ν = L ν /4πD 2 , where D is the distance to Earth.Thus, 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) 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.

Appendix B. Equipartition
The equipartition calculation follows.As shown above, jet kinetic power is where dm dt = ρ dV dt = ρA dx dt = ρAu Kinetic energy density is [6]  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 (Figures A2 and 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.

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 climbcurrent, representing the current LOS advance in the x ′ , y ′ , and z ′ axes, respectively (Figures 1 and 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).We then proceed to calculate the current length of the LOS where LOS length is measured in cell length units and nx1current = nx10 + rc, ny1current = ny10 + uc, nz1current = nz10 + cc (A16) 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).where Γ π is the pion Lorentz factor and t π0 = 2.6 × 10 −8 s (A33)

Appendix C.3. Relativistic Effects-Doppler Boosting
. 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 π .
In the calculations, we employ the form 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
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: where⃗ r is the spatial position vector, N is the particle density of the daughter population, and Q is its injection function.

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.
x is the fraction of the pion energy to proton energy, and n(z) is the jet flow proton density.versus pion energy E π .In order to obtain the pion distribution, the following transport equation is employed: where N π (E, z) denotes the pion energy distribution.Then and ) 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 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.
1 c (f pp(4d) = 1.0 if |(u/u th )| > 1, and f pp(4d) = |(u/u th )| z fpp if |(u/u th )| < 1); f obs is the observing frequency; and f obs D 4D local is the calculation frequency (blue-shifted or red-shifted according to the relative orientation between the local LOS and local velocity); − → B 4d is the local magnetic field vector; − − → LOS is the line of sight (LOS) vector at an angle ( − − → LOS, − → B 4d ) to the magnetic field spatial vector; D 4D local is the local Doppler factor of E/M emission from a computational cell, calculated as D = √ 1−u 2

Figure 2 .
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 (Table1).

Figure 3 .
Figure3.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 (Table1).

Figure 4 .
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.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 6 .
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 8 .
Figure8.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 9 .
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 11 .
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 12 .
Figure12.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 Figure10.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 A2 .
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 A3 .
Figure A3.Simultaneous advance in two-dimensional space and forward in time of a few lines-ofsight.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 A5 .
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 A6 .
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.
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 −1 loss = t −1 sync + t −1 adb + t −1 pp (A36)The time-scales of the different energy-loss mechanisms are presented in FigureA8.

Figure A9 .
Figure A9.Density of nonthermal protons in the jet using a high-energy cutoff feature plotted with energy.

Figure
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 πEp .Calculated at three different energies of the hot proton.

Figure A11 .
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 A13 .
Figure A13.Flow diagram of nemiss depicting the basic structure of the programme.

Figure A14 .
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 A15 .
Figure A15.Flow diagram of rlos.1 depicting procedures and functions called during programme execution.

Table 1 .
Three different imaging runs based on three separate underlying hydrocode runs.