Synthetic Neutrino Imaging of a Microquasar

Microquasar binary stellar systems emit electromagnetic radiation and high-energy particles over a broad energy spectrum. However, they are so far away that it is hard to observe their details. A simulation offers the link between relatively scarce observational data and the rich theoretical background. In this work, high-energy particle emission from simulated twin microquasar jets is calculated in a unified manner. From the cascade of emission within an element of jet matter to the dynamic and radiative whole jet model, the series of physical processes involved are integrated together. A programme suite assembled around model data produces synthetic images and spectra directly comparable to potential observations by contemporary arrays. The model is capable of describing a multitude of system geometries, incorporating increasing levels of realism depending on need and available computational resources. As an application, the modelling process is applied to a typical microquasar, which is synthetically observed from different angles using various imaging geometries. Furthermore, the resulting intensities are comparable to the sensitivity of existing detectors. The combined background emission from a potential distribution of microquasars is also modelled.


Introduction
Microquasars (MQ) comprise a binary stellar system where a main sequence star orbits a compact object, either a neutron star or a black hole [40].Matter from the star accretes onto the collapsed stellar remnant, resulting in the production of twin relativistic jets pointing in opposite directions.Those jets emit over a broad spectrum, from radio to very high-energy (VHE) γ rays and neutrinos [41][42][43][44][45][46][47][48].
As mentioned in [41], apparent superluminal motion in certain MQs indicates the presence of bulk hadron flows in the jets.The assumption of equipartition [45] leads to high magnetic field estimates for the jet [49].This, coupled with the fluid approximation for the jet matter due to the presence of tangled magnetic fields [50,51], allows for magnetohydrodynamic (MHD) approximation for the jets.A toroidal magnetic field component may retain jet collimation over considerable distances along its path [49,52].Moreover, external confinement from surrounding winds is equally important [45,53].
In order to study the jets, a selection from among the wealth of theoretical results is compared to observations of those remote systems.The relative scarcity of detailed data is complemented by the use of numerical simulations of a jet system, where a model setup is evolved and then imaged synthetically.The final model emissions are placed next to observations, running many examples until a match is achieved.If no positive detections exist yet, then a general match to theoretical results and the sensitivity of active observing arrays is desired.As a next step, going backwards, the jet model is reverse-engineered to its initial boundary and generally internal or unobservable conditions, which emerge as the link between jet theory and observations.The above process can offer increased insight into the inner physical workings of the jets and their surroundings, allowing for their study as a complex, evolving dynamical system.A more accurate description of the system of interest is then obtained.
In this paper, the production of VHE neutrinos from generic MQ jets is modelled using the method of dynamic and radiative relativistic MHD simulation.A set of surrounding winds assists with the confinement of the jets, adding realism to the model.Within the jets, a complex turbulent environment allows for the production of a variety of different signals, from radio to X and γ rays.Furthermore, cascades of highenergy particles produced in the jets lead to an ecosystem of different particle populations connected through transport phenomena.The production of neutrinos that leave the system opens the possibility of detection on Earth from modern arrays.
The solution of the transport equation from one particle distribution to the next, along a cascade, allows for the expression of the intensity of emitted neutrinos as a function of dynamic and radiative jet parameters at a given point.This way, local model parameters at each space-time point in the model jet are directly connected to the final particle emission at the same point.Repeating the latter process for a number of energies provides a neutrino energy spectrum at each jet space-time point.Line-of-sight integration follows, leading to the production of a synthetic neutrino image of the system and a whole-jet neutrino energy spectrum.
The paper is organized as follows.In Section 2, the theoretical background of the work is presented.In Section 3, the emission of particles from the jet is obtained.In Section 4, results are presented and discussed.Normalization and equipartition (and the synthetic imaging process) are described in Appendixes A and B respectively.

Theoretical Setup
In our generic MQ model, an accretion disk is assumed around the compact object [54].Twin jets emanate from the vicinity of the collapsed star, collimated by a toroidal magnetic field component.Adopting a heavier pair of jets, their kinetic power was set to L k = 2 × 10 38 ergs −1 (see Appendix A).The authors of [45] argued a 10% Eddington luminosity jet power, leading to L k = 10 38 ergs −1 for a 10 M ⊙ black hole, which is comparable to our simulation.Furthermore, for the ratio of proton jet power L p to electron jet power L e , the same authors argue either L p L e ≃ 100 or ≃1; we adopted the former hypothesis, favouring protons.As a first implementation, we calculate neutrino emission originating from the influence of the high-energy proton distribution, while there is also a potential comparable contribution from the corresponding high-energy electron distribution [45].
In the jets, equipartition is assumed between kinetic (ρ k ) and magnetic (ρ B ) energy densities, meaning ρ k = ρ B ; therefore, at each jet point ⃗ rz, the CGS magnetic field is B(⃗ rz) = 8πρ ⃗ rz [50,51], a close match with the B used in the simulation (see Appendix B).External magnetic fields tend to be quite smaller [55]; therefore, as a first-order approximation, they are not included in the surrounding winds.

Non-Thermal Proton Density
Neutrino emission from the jets is taken to originate from proton-proton interaction between a distribution of hot (fast) protons and cold (bulk flow) protons [41,44,45,48,56,57].Some of the bulk protons are accelerated at shock fronts according to the first-order Fermi acceleration mechanism, with a time scale of [50,51,58] where B is the magnetic field, E p is the proton energy, e the proton charge, and c the speed of light.η = 0.1 represents an acceleration efficiency parameter, assuming efficient acceleration in moderately relativistic shocks in the vicinity of the jet base [58].
As an approximation, high-energy electron distribution is deferred to future work.Focusing on hadrons, we adopt a power-law distribution for the relativistic protons as a function of their energy E of the form N p = KN p(0) E −α [53], K being a scaling constant connecting fast hot and thermal proton densities, N p being the bulk proton density at a given jet point and N p(0) the density at a reference jet point, with either α ≈ 2 [44], or a variable α [48], where α is the proton spectral index in the local jet cell matter frame.Alternatively, a transport equation could be used to find the distribution [44].
As a further approximation, the aforementioned hot proton distribution is taken to be isotropic in the jet frame, assuming that, at each jet point, l sc < l r , where l sc is the scattering length and l r the radiative length, a hypothesis backed by the nature of diffuse shock acceleration [59].

A Note on Jet Frame Anisotropy
For protons accelerated at diffuse shocks, the above assumption of isotropy is justified by the need to preserve, after every bounce, at least some proton energy [51].Consequently, scattering length l sc is less than radiative length l r .Otherwise, the proton would not have any energy left after the bounce, negating the acceleration process.
According to the work in [59], assumed anisotropy of hot proton distribution can be reflected to neutrino distribution.In the jet system, emission would then be projected off axis, even reinforced, under certain orientations, in the lab frame.

Proton Energy Loss
Following the works in [45,48,56], certain energy loss mechanisms are included.This presentation is performed in a cell with the properties of (u x , u y , u z , b x , b y , b z , n, ϕ 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).In the latter, u stands for velocity and b for magnetic field along directions x, y, or z.Bulk flow proton density is denoted by n, and ϕ1 ≃ ϕ (ϕ 2 = 0) is the complementary to θ, the angle to the line of sight.Last, α is the high-energy proton distribution spectral index.As an exception, the pion injection function presentation uses a different velocity vector of (0.2, 0.8, 0.1)c.Nevertheless, in the model runs, these are potentially performed in every cell.
We consider cut-off E for protons E ≤ 10 6 GeV.For the adiabatic expansion time scale, we have [45] where z j = 10 11 cm is the characteristic lateral size scale of the jet.For this simple calculation, u b(adb) is preset to 0.8c.For the p-p collision loss mechanism, we have where n is the bulk flow proton number density, K pp = 0.5 [45] and σ is the inelastic p-p collision cross section [45] σ where E th = 1.2 GeV and L = ln(E p /1000 GeV)) (see in [45,56]).Equation (3) is justified if we consider a small cube of matter of number density n, moving at speed (near) c and 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 σ  For the pion decay time t π and the characteristic pion decay timescale t π0 , we have the following equations: t π0 = 2.6 × 10 −8 s (5) and where Γ π is the pion Lorentz factor, which, in practice, takes the form (m π is the pion mass and E π the pion energy) where light escape time t esc strongly affects the final result.The synchrotron loss time scale is defined by [45] t 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, e is the electron charge and B is the local magnetic field.The form of the latter term E p m p c 2 , which is equal to the proton Lorentz factor, which is essentially Γ p , facilitates energy-dependent calculations later.In total, In Figure 2, the various energy-loss mechanism time scales are presented.

Model for the Interaction of Thermal and Non-Thermal Protons in the Jet
Hot-cold proton interaction results to a distribution of high-energy pions, which then decay, allowing for the creation of energetic neutrinos.We have [60][61][62][63] pp → ppπ 0 + π 0 , ( for neutral pions π 0 , and for π ± .π 0 decay to gamma rays, while π ± mostly decay to an antimuon or muon and a muonic neutrino or antineutrino (prompt neutrinos) [60,63].
As an approximation, we neglect both neutrino production through secondary channels and delayed neutrinos.
For each successive particle population in the above cascades, the transport equation [64] can be solved.
The transport equation for nonstochastic phenomena and for time-independent transport (transport time much less than the time step of the dynamic simulation) takes the following simplified form: ∂N ∂E where t loss is the decay timescale for the particle in question, ⃗ r is the location vector in space, N is the particle density of the produced particle population, and Q is its injection function.Q is calculated from the previous population up the cascade.For example, N may represent protons and Q, (which includes proton acceleration effects), is expressed as a function of the local thermal proton distribution.In turn, the N of pions can be obtained using their Q, which in turn is a function of the hot proton N, and so on along the cascade.Nevertheless, a power-law distribution is assumed for protons, skipping having to solve the first transport equation in the cascade.From protons to pions, then to muons and neutrinos, each generation of particles leads to the next one.The authors of [56] calculated the properties of resulting particle distributions over a large energy range, performing Monte Carlo calculations with the results of particle physics.

PLUTO Code
PLUTO [65] 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 carried out 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 programme may run efficiently in parallel on various platforms.
In previous works [62,66,67], the hadronic jet was modelled using the PLUTO code.PLUTO results were then processed in order to calculate the emissivity of γ rays and neutrinos using various approximations.As far as neutrinos are concerned, emission was calculated at only a handful of points along the jet, thus not taking advantage of the detail offered by a numerical jet simulation.In this paper, the emissivity of neutrinos is separately calculated at each computational cell using the angle (los,u) formed between LOS and local velocity.This calculation is huge compared to the previous one, but the benefit is that a result is obtained at each point.New code NEMISS [68] (not available before) is employed here, which performs the calculation on data produced by PLUTO.
Furthermore, a time-delay-capable line-of-sight code, RLOS2 [69], is employed here, which was not available in the aforementioned previous works (much simpler LOS code [70] used back then was also written by this author).RLOS2 reads the combined results of PLUTO and of NEMISS, and produces synthetic neutrino images of the model system using either a focused beam geometry or parallel rays.
The improved relativistic transformation of the hot proton distribution by [71] is now employed as opposed to [72] in the previous works.
The PLUTO jet is now a twin in 3D space, using new files for setting up PLUTO for our problem [73], adding to the realism of the system.The twin jet model system employed here is an evolution of the single jet system, also set up by this author in PLUTO [74].The magnetic field is now adjusted for equipartition, and the model parameters are generally more refined compared to those in earlier works.

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) is found in [71,72].
The latest variant originates from [71] (TR), used here, minus a geometry factor that we absorb into the normalization factor where Γ is the Lorentz factor of the particles in a tiny volume.The particle mass is m and β = u/c is the particle speed in units of the speed of light c.E is the particle energy and α is the spectral index of the distribution.
On the other hand, the authors of [72] (PS) say Equations ( 14) and ( 15) are compared in Figures 3 and 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 π [45,48,56].
where x = E/E p . Figure 5 shows p-law fast proton density.Figure 6 xF π 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 x of all potential p-p interactions.
x is the fraction of the pion energy to proton energy, and n(⃗ rz) is the jet flow proton density.In order to obtain pion distribution, we solve the following transport equation: where N π (E,⃗ r) denotes the pion energy distribution.We proceed where The quantity τ π (E ′ , E) is the pion optical depth and b π is the energy loss rate of the pion.As an approximation, the last term in the latter expression is omitted.Figure 8 plots NEMISS software function U (U analytical ), representing N π , with pion energy.The above are performed for each computational cell, where the quantities for radiative purposes are considered locally constant.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 the radiative interactions belong to a single time step each time.

Neutrino Emissivity
The emissivity of prompt neutrinos [44,45,56,57] is where E is neutrino energy, r π = (m µ /m π ) 2 , x = E/E π , and t π is the pion decay timescale.Θ(χ) is the theta function [45,62].Neutrino emissivity is calculated for each individual cell using the cell's own angle to the LOS crossing that 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 [69].A synthetic image of the model system is thus produced.

Results and Discussion
Using the formalism presented in this paper so far, the neutrino emission at each computational cell of the model is calculated.This method is heavier from a computa-tional point of view, but allows for obtaining a separate neutrino emission from each spatiotemporal point of the twin jet model.Thus, we aim for the result where intensity I is calculated at the 3D computational cell at⃗ r, represented by the x, y, and z coordinates of the cell.Time t is obtained from the time tag of the PLUTO data dump where the cell belongs.Thus, the above equation is globally applied to all selected PLUTO data (the user may select beginning and end times for the global calculation).We now proceed to describe the setup of the simulation.The jet base is situated near the centre of a Cartesian grid.A continuous model jet representing a microquasar system is injected at a u jet = 0.865c (a Lorentz factor of 2, which lies between higher microquasar Lorentz factors used in the literature, such as Γ = 5 in [41], and Γ = 5/3, corresponding to u = 0.8c, a characteristic value for the jets in GRS1915 + 105) is studied with the RMHD setup of the PLUTO hydrocode, at a uniform grid resolution of 60 × 100 × 50.Grid size is (120 ×10 10 cm) × (200 ×10 10 cm) × (100 ×10 10 cm); therefore, cell length is 2 ×10 10 cm.The grid size is such that it focuses on the area of the inner jet, where γ ray and neutrino production is expected.This minimal size of the cell means that a starting radius of the jet of a few times 10 10 cm is necessarily implied, as a few cells' diameter of the nozzle is used.This compromise is imposed by the nature of the employed simulation, which utilises a homogeneous grid.In future work, a non-homogeneous grid may allow for better focusing on more realistically resolving the jet input nozzle.
In all of the model runs, the same initial jet density of 10 10 protons/cm 3 was used (a typical value for the inner microquasar jet, also compatible with the energetics of the jet and its kinetic luminosity), 2000 times less than the maximal surrounding gas density (i.e., a light jet is assumed, which is a possibility that supports a rich jet-wind interaction environment, but is also more demanding from a computational point of view).Winds comprise an accretion disk wind construct and a stellar wind that falls off away from the companion star, located off-grid at (4 × 10 12 cm, 1 × 10 12 cm, 4 × 10 12 cm), while the jet is threaded by a strong confining toroidal magnetic field of B = 10 4 G, assuming equipartition between kinetic and magnetic energy density (see Appendix B for the calculation of the latter equipartition value for B, in relation to the jet kinetic luminosity).A guide for inner system winds and their densities was SS433 [54].Simulations were run until t = 842 s, saving a data snapshot every 25 (simulation) s.A three-dimensional snapshot of density is shown in Figure 9, where we can see the magnetically collimated jet pair advancing through surrounding winds.
The above figures show a narrow jet barely expanding into its surrounding winds.This small half-angle is then rather counterintuitively expected to result in a faster decline of neutrino emission with energy, as discussed in the discussion section of [45].
A number of empty user parameters were employed in order to house particle emission results later.Then, the above PLUTO run was copied into many directories.In each, the NEMISS programme [68] was run, which calculates neutrino emissions for a specific imaging geometry and setup.This programme is able to read 4D spatiotemporal data output from PLUTO into a 5D array, which also includes particle energy as a fifth dimension.Then, NEMISS calculates the neutrino emission at each point of the 5D data array.Results were overwritten into suitably prepared data files of the originally empty user parameters of the hydrocode.Thus, NEMISS processes PLUTO output to include a neutrino emission spectrum at each spatiotemporal data point.PLUTO data processed by NEMISS are then ready to be read by relativistic time-delay LOS imaging programme RLOS2 [69], which produces synthetic neutrino images of the system.Over a string of particle energies, the intensity sum of the whole synthetic image of the jets is calculated for each energy.Thus, the plot of jet neutrino intensities is produced.
The intensity plots of the model pair of jets are created using Veusz, a software for plotting data written by Jeremy Sanders and contributors, and distributed under the GNU/GPL licence.RLOS2 and NEMISS were written by the author and are available under the lGPL licence.PLUTO was written by Andrea Mignone and collaborators, and is available under GNU/GPL.
Table 1 shows a number of simulation parameters.Those include computational cell length, jet density, and both winds' maximal densities (those gradually declined away from their sources).In PLUTO, the piecewise linear method was set up using the MUSCL Hanckock integrator.An ideal equation of state was used.The binary companion is located outside the grid, and was estimated to be at most up to an order of magnitude greater than that of the compact object.Jet speed is 0.866c, while its kinetic luminosity is 2.5 × 10 38 erg/s.As a first use of the programme suite, a rather low spatial resolution of 60 × 100 × 50 was employed in PLUTO in order to accommodate for the heavier neutrino emission calculation later.As far as RLOS2 is concerned (synthetic imaging), either focused beam or parallel rays are employed as an imaging method, while the time-delay effect of RLOS2 is not employed at this stage, as it requires multiple RAM memory to be used properly.The synthetic image is projected either on a side of the computational box, either front or sidereal, or on a fiducial imaging screen, again either frontal or sidereal.More specifically, a series of imaging geometries were employed following the RLOS2 programme convention: (imaging geometry) Case 1, parallel rays projected onto the XZ plane; Case 2, the same but onto the YZ plane; Case 3, focused rays onto the XZ plane; and Case 4, focused rays onto the YZ plane (see also Table 1).Three different angles were employed for Cases 1 and 2, while for Cases 3 and 4, respective focal points implied near-head-on and sidereal views.
RLOS2 [69] was then run using the combined PLUTO-NEMISS data with sfactor = 1 for the pload shrink factor.In general, the imaging process may or may not use all snapshots available to it depending on the light crossing time of its model segment (adjusted through the clight parameter in RLOS2).Trying to read more snapshots than what is loaded corrupts the hydrocode time array of RLOS2, called T, resulting in errors.For simplicity, in our case, an artificially very high clight was used in order to effectively switch off the time-delay effect.A double filter was used for velocity and for los,u angle.A minimal velocity and maximal angle were set in order to trigger the calculation of the neutrino emission for a particular cell.This way, the expensive part of the simulation was only performed where it was really worth it.This partly alleviated the discrepancy between computational costs of the dynamic and the radiative parts of the model.
The twin jet simulation used in this work represents a single fiducial microquasar using characteristic properties.This system was dynamically set up to be relatively close to a number of microquasars, such as Cyg X-1 or GRS1915 + 105.From this point on, the model system is imaged with different methods and at different angles in order to explore the perhaps dominant effects of orientation, both locally and globally in the jet.Those imaging results can then be extrapolated to a variety of similar microquasar systems.
An important aspect of this modelling approach is that each cell has different visible emissivity from Earth than that of its neighbours.That is because each cell may differ from the next one in terms of both speed and orientation to us.This combination means that the hydromodel generally gives different results than those of the steady-state one.A vortex with relativistic velocities, for example, may partly appear very luminous where it is fast with local speed pointing towards us, and also too dark where velocities point away from us.In this simulation, such effects were limited, but at a higher resolution, it is expected that nonlinear dynamic effects in the hydrocode profoundly interact with the radiative part of the model.
The scale of the total emission increases the closer that the LOS approaches to the jet pair axis (Figure 10).The employed low resolution did not allow for significant nonlinear dynamic effects to appear, yet the concept of the modelling process was proven to work in its entirety.On the other hand, the normalization process demonstrates the possibility of potential observations, as the results potentially fall within the detection range of contemporary arrays [45].The detection ability of km 3 array is depicted in the normalized spectral emission distribution (SED) plots, as a measure of comparison with the model results.A certain potential for detection appears that is rather promising to explore.
Figure 10.Normalized SED from a series of radiative simulations where comparisons with potential observations are possible.See Appendix A for more details on the normalization process.Angle ϕ 1 is complementary to the angle between jet and the LOS, θ (ϕ 1 = 90 − θ).Then, angle ϕ 1 is nearly 90 degrees when looking along the jet axis.Therefore, ϕ 1 ≃ 90 degrees in imaging geometry's Case 3 (rays near parallel to the jet axis), and ϕ 1 ≃ 0 in imaging geometry's Case 4 (rays nearly perpendicular to the jet axis).Intensity decreases as the angle of observation moves away from the jet pair axis.For comparison, the KM3NeT threshold [45] for detection is included as a dotted line near the top.Image produced with Veusz using data produced with PLUTO, RLOS2, and NEMISS.
More specifically, we can see in Appendix A that intensity on Earth is proportional to kinetic jet luminosity L k and inversely proportional to the square of the distance to us D 2 .Consequently, a sample set of rates can be extracted from the model and used as a reference for other microquasars at different distances and with different jet energies than those of the standard.Figure 11 shows the weighted set of rates expected on Earth for a sample microquasar viewed at 30 degrees from the jet axis with L k = 10 38 ergs −1 and D = 5 kpc.Other systems then have , where I, L k , and D refer to a new microquasar; I 0 , L k0 , and D 0 represent the standard plotted here, and the profound effect of the viewing angle is implicitly included.The result of Figure 11 is comparable to the sensitivity of the state-of-the-art instrument arrays IceCube/KM3NeT IceCube in terms of a squared energy weighted curve, which falls below 10 −8 throughout the plot's energy span [75].
The above estimate may then be employed in order to provide a rough estimate of expected neutrino emission from a distribution of microquasars in the galaxy.The authors of [76] argued an estimated population of approximately one-hundred systems in our galaxy.Furthermore, their discussion of γ ray emission from microquasars clarifies 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.
We proceed by accepting 100 systems at various distances ranging from a minimum of 1 kpc to a maximum of 30 kpc, with average kinetic luminosity similar to our model system.The linear dependence of emissions 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 a system situated at 30 kpc has 36 times less than that.Last, an orientation of less than 60 degrees might be 1 order of magnitude less than our value, but a jet system aimed towards us could have up to 100 times more visibility at Earth unless a very fast jet occurred.Consequently, the single most important factor is orientation, followed by distance and lastly by jet kinetic power.The latter order allows for an estimate of perhaps 5%, or five systems with a very high relativistic boosting towards us, a number of maybe 40 or 50 at angles above 45 degrees, and lastly maybe 50 at below 45 degrees.The first five probably contribute the most on average, and the ones viewed from the side have a smaller effect.A possible system at a smaller distance would of course dominate the distribution, but the possibility for such an occurrence is questionable.
On the basis of the above discussion, we then accept 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 30 degrees from the line of sight, which is the case used in Figure 11.The reason for having the average angle at less than 45 degrees is the higher contribution from systems aimed towards as.We then multiply our single microquasar result by 100 (population size), divide it by 3 2 (distance) and leave the jet power effect at unity.A rough first estimate could then be to multiply our single system result at 30 degrees from the jet axis by a factor of ten (Figure 12) and then use it for comparison with observations.
Orientation seems to play a crucial role here and is thus given the primary role in the synthetic imaging process by employing various orientation scenarios for the model pair of jets.In addition, this model calculates the effects of orientation at each point of the 3D PLUTO twin jet simulation, and then produces a synthetic neutrino image.Thus, the important effects of differential projection effects are explored, paving the road for more detailed simulations in the future using this programme suite.
In contrast, previous similar works [61,62], calculated neutrino emission at just a handful of points along a single model jet (a much smaller computational task), and then used a semi-analytic approach to cover the rest.Furthermore, a number of programme improvements were incorporated into the models, such as using [71] for relativistic orientation and velocity transform of the hot proton distribution, as opposed to [72] in the previous works where this author contributed.
The above results for microquasar distribution may vary to either direction by possibly an order of magnitude, subject to a more detailed statistical analysis.This is because there are similar systems with higher or lower jet kinetic power, as well as systems with various individual properties.Nevertheless, it seems possible that the detection of a background emission from a potential distribution of microquasars in the galaxy is within the realm of modern detector arrays.This is also a consideration for the next generation of new or upgraded arrays being planned today.On the other hand, a single X-ray binary system also looks promising as a galactic source of high-energy neutrinos.This is a potential target for a particle sensor with 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 working in electromagnetic spectra, such as radio, X-rays, and γ rays.In such a case, a neutrino observation of a microquasar may form part of a multi-wavelength observation effort aimed at the system of interest.Intensity on Earth at an angle ϕ 1 = 60 degrees, meaning around 30 degrees from the jet axis, weighted by the energy/energy squared of the particle.These represent bin plots, as each data point lies on a higher size scale than the next one.Furthermore, they are also cumulative rate plots, upwards from each given energy, as rate contributions from higher energies are much smaller than the starter one.In comparison with the sensitivity of IceCube/KM3NeT of below 10 −8 in this plot, results seem marginally acceptable in anticipation of potential detection.Image produced with Veusz.For comparison, the KM3NeT threshold for detection is included as a dotted line near the top.

Final Remarks
Particle emission from a typical microquasar was simulated using a suitable programme suite.Results verified the integrity of the process, paving the way for more detailed runs.Furthermore, the model was employed in order to provide particle emission estimates for both a single microquasar and a potential galactic distribution of such systems.The latter approach facilitates a comparison with the output of contemporary detection arrays, where microquasars could contribute to a background of high-energy neutrinos.
In the model, a series of both dynamical and imaging parameters may be adjusted in order to cover different scenarios.The programme suite works in a highly automated manner, and is prepared to take on higher-resolution applications where the relativistic effects of nonlinear dynamics may appear in full.
The ability to focus on individual cells could greatly differentiate each jet element from the next in terms of emission.An MHD jet has great local variability in both particle and radiation emission intensity in any given direction.The detailed dynamics of the jet influence the appearance of the system depending on both the direction and magnitude of the local velocity, and on pressure and density.Consequently, a jet system with turbulence, vortices, colliding with clouds, etc. is expected to be subject to the aforementioned local variations of intensity.
As far as absorption is concerned, the model may directly include the emission and absorption of electromagnetic radiation at different frequencies.Should adequate computing resources be employed, the time-delayed description in the programme can also be activated.For example, a turbulent relativistic jet colliding with a cloud has different parts of it moving at high velocities in different directions.The image is then dynamically formed, the rays crossing a choreography of relativistically moving jet elements.The final image may be quite different than what is initially expected, as demonstrated, for example, by the effect of apparent superluminal motion.For a complex jet system, running the model at higher resolutions with the time-delay module could reveal many physical details, drawing a more realistic picture of the system.
In general, microquasars may locally emit at reinforced levels of intensity due to the combination of jet dynamics and relativistic projection.The reason can be internal jet turbulence or interaction with clouds and surrounding winds.For γ rays and neutrinos, such dynamic effects should occur in the vicinity of the jet base.
Furthermore, the employed model can be used as a basis for expanding the approach to systems of different scale.The innermost AGN jets can be sources of ultra high-energy cosmic rays [77], and PLUTO can model those jets with a suitable set of initialization parameters.Special relativistic MHD should be employed as an approximation, though.There is a possibility to include a quasi-Newtonian potential as an improved approximation for the innermost part of a quasar jet.The emission model, which in our case was NEMISS, should be altered in order to include the new emission physics.Synthetic imaging code RLOS2 is ready to use with any emission and absorption input, and only minor changes are required.
Further out along an AGN jet, neutrino emission may occur from high-energy proton acceleration along with other signals such as γ rays [78].This description is similar to microquasars, and only the scales differ.Consequently, it should be possible to suitably adapt the current simulations in order to model neutrino emission from the inner part of a quasar jet.Data Availability Statement: Synthetic data used in this paper were produced with code linked to in the bibliography.cubical computational cell at the jet base.We then express density as a function of proton number density N p and proton mass m p , ρ = N p m p . (A6) 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 L k ), L ν = αL k , representing the portion of jet power emitted in neutrinos.For normalisation, a working value is taken as 10 −3 .This can be justified from a q rel = 0.1 for the energy content of the relativistic particles in the jet [45,48], on top of which we employ the efficiency of the cascade when transferring energy from hot protons to final neutrinos.
The shape of the spectrum is also affected by acceleration efficiency [48], and from the opening angle of the jet [45], 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 non-thermal protons to the neutrinos.
Furthermore, we introduce a factor α = L ν /L k , representing the portion of jet power emitted in neutrinos.A typical value is taken as 10 −3 .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 The intensity of the jet is then expressed as I ν = L ν /4πD 2 , where D is the distance to Earth.Thus, In our simulation, the jet beam travels at β = u c =0.866, with a density of 10 10 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.Furthermore, Γ = 2. 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 an arbitrary units neutrino intensity plot, for the case of nearly non-beamed data, at ϕ 1 = 10 degrees.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, admitting 10% coverage per order of magnitude scale level.Thus, we find about 10 11 , which means that our sum is 10 times smaller, or approximately 10 10 , in (AU)*GeV, where AU stands for arbitrary units.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 unnormalized intensity plot with neutrino energy, expressed in units of C 0 , in order to find the latter (normalization constant) (Γ(N p m p )N cell L 2 cell )β 3 c 3 = (PLOTAREA) * C 0 erg/(s cm 2 ) GeV (A9) For our case, we find C 0 ≃ 2 ×10 −21 , 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 10, which may be directly compared to other models and to observations.

Appendix B. Equipartition Calculation
The equipartition calculation now follows.As shown above, the jet kinetic power is where dm dt = ρ dV dt = ρA dx dt = ρAu

(
inel pp ) pp , resulting in the inverse time scale for the aforementioned p-p collision.In Figure 1, σ (inel) pp is plotted.

Figure 1 .
Figure 1.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 2 .
Figure2.High-energy proton distribution loss time scales, for various processes in the jets, plotted 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 3 .Figure 4 .
Figure3.Ratio of high-energy proton distribution density transformation as calculated by the formulae of TR and PS, respectively, for three different angles.TR stands for[71], PS for[72].A reproduction, for verification, of a figure from in[71].

Figure 5 .
Figure 5. Density of non-thermal protons in the jet using a high-energy cut-off feature plotted with energy.

Figure
Figure 6.F (pp) π x, E x function, Equation (17), F function, corresponding to the pion spectrum emerging from a single (hot-cold) proton collision, multiplied by the x = E π Ep fraction.Calculation performed at three different energies for the non-thermal proton.

Figure 7 .
Figure 7. Pion injection function Q, weighted by pion energy, measured in non-normalized units, describing the combined spectrum from a multitude of (hot-cold) p-p collisions.We can see contributions rapidly declining as particle energy increases.As an exception, this figure uses a velocity vector of (0.2, 0.8, 0.1)c.

Figure 8 .
Figure 8. Pion energy distribution plotted in non-normalized units versus energy.In the software, the above distribution is represented by function U.

Figure 9 .
Figure 9. Three-dimensional side view of twin model jet system.Snapshot 14 of u = 0.866c hydrocode run corresponding to a model time of t = 350 s (14 × 25), depicting the density in a logarithmic plot.Both jet fronts are advancing towards the ends of the grid, traversing the surrounding stellar wind after crossing the simplified accretion disk wind construct.Image produced with VisIt.

Figure 11 .
Figure 11.Intensity on Earth at an angle ϕ 1 = 60 degrees, meaning around 30 degrees from the jet axis, weighted by the energy/energy squared of the particle.These represent bin plots, as each data point lies on a higher size scale than the next one.Furthermore, they are also cumulative rate plots, upwards from each given energy, as rate contributions from higher energies are much smaller than the starter one.In comparison with the sensitivity of IceCube/KM3NeT of below 10 −8 in this plot, results seem marginally acceptable in anticipation of potential detection.Image produced with Veusz.For comparison, the KM3NeT threshold for detection is included as a dotted line near the top.

Figure 12 .
Figure12.Intensity on Earth, of a fiducial distribution of systems weighted by the energy/energy squared of the particle.Again, these may represent bin plots due to the logarithmic decline of the quantity.For the same reason they are also cumulative rate plots, upwards from each given energy.For comparison, the KM3NeT threshold for detection is included as a dotted line near the top.Image produced with Veusz.

Funding:
This research received no external funding.Institutional Review Board Statement: Not applicable.Informed Consent Statement: Not applicable.

Table 1 .
Five different imaging runs based on same underlying hydrocode run.