Neutrino Emission from Jet Interactions in Microquasars

1. Introduction

Microquasars (MQs) are binary stellar systems comprising a main sequence star orbiting a collapsed stellar remnant [1]. Mass transfer onto the compact object powers relativistic jets launched largely perpendicular to the orbital plane. These jets emit radiation over the electromagnetic spectrum, from radio to very-high-energy (VHE) gamma rays, and also neutrinos [2,3,4,5]. Observations such as apparent superluminal motion support bulk hadronic flows in these jets [2].

Recent ultrahigh-energy (UHE) gamma-ray observations strengthen the case for Galactic compact-object jet systems as potential PeVatron accelerators. In particular, results from LHAASO provide a new UHE view of Galactic gamma-ray sources and motivate renewed interest in jet-powered systems as extreme particle accelerators [6]. Furthermore, UHE gamma-ray emission associated with black-hole jet systems has been discussed in the context of LHAASO detections, providing additional motivation for hadronic acceleration scenarios and associated neutrino production [7]. In parallel, microquasar jet and jet–environment interaction models have recently been explored as PeVatron candidates, including jet–cocoon systems and microquasar remnants as potentially “hidden” PeVatrons [8,9,10].

These developments motivate time-dependent, simulation-based modeling that links RMHD jet dynamics and interaction sites (shocks, plasmoid fronts, terminal regions) to neutrino observables. Since MQ distances and complex geometries limit direct inference, detailed numerical modeling is essential to connect theory with multi-messenger observables.

Strong magnetic fields near jet bases and their tangled structures justify a fluid description via special relativistic magneto-hydrodynamics (RMHD) [11,12,13,14]. Toroidal magnetic components collimate the jet [11,14], while interactions with stellar and disk winds shape jet morphology and confinement [4,15].

Our study combines RMHD jet evolution with neutrino emission modeling, including jet–environment interactions and relativistic effects (beaming, time delays). The focus is on neutrino emission enhancement due to plasmoid collisions with ambient photon and matter fields, with explicit statements of acceleration zones, diffusion regimes, magnetic-field prescription, target photon fields, and the relative roles of $pp$ and $p\gamma$ processes.

The paper structure is: Section 2 covers particle acceleration and radiation physics (including explicit zone prescriptions), Section 3 details neutrino emission formalism and target fields, Section 4 describes software tools used, Section 5 provides numerical setup and parameters, Section 6 presents synthetic neutrino results and observability including distance-dependent detectability, and Section 7 summarizes main findings.

2. Theoretical Background

2.1. Particle Acceleration Regions, Mechanisms, and Diffusion Regimes

High-energy protons are accelerated at shock-like structures that naturally arise in intermittent relativistic jets. In the present work we adopt a zone-based description to explicitly specify where acceleration is assumed to occur and what transport regime is used.

Zone A: Internal shocks / plasmoid fronts. Internal shocks form due to velocity irregularities and intermittency; in addition, plasmoid fronts and compressed regions appear due to jet instabilities and jet–wind interactions [3,13]. In the RMHD output, candidate acceleration cells are identified by compression and shock proxies (e.g., negative velocity divergence $\nabla ·\mathbf{v}<0$ together with elevated pressure gradients). These regions typically dominate the time-dependent nonthermal power injection.

Zone B: Recollimation / interaction layers. As the jet propagates through ambient stellar and disk winds, recollimation and shear layers can form. These regions can host additional shocks and turbulence, contributing to acceleration and re-acceleration.

Zone C: Terminal / matter-loaded regions (head/cocoon). Toward the terminal regions of the jet where matter accumulates, thermal densities can become high, and hadronic $pp$ interactions may become significant. This is also where transport may differ from internal plasmoid fronts due to evolving turbulence levels and larger coherence scales.

Acceleration mechanism. We assume first-order Fermi shock acceleration in the above shock-like regions. The acceleration timescale is modeled as

$$t_{\mathrm{acc}}^{-1}=\eta {\displaystyle \frac{ceB}{E_p}},$$

where $\eta$ is an efficiency parameter (order $0.01$—$0.1$ for relativistic shocks), B is the local magnetic field, and $E_p$ the proton energy.

Diffusion regime and region-dependent transport. Particle transport is controlled by the turbulence level and magnetic-field structure. We adopt an explicit diffusion coefficient parameterization

$$D(E)=D_0{\left({\displaystyle \frac{E}{E_0}}\right)}^{\delta },$$

with $\delta$ chosen by zone to reflect different turbulence regimes. In strongly turbulent, shock-compressed plasmoid fronts (Zone A), we use near-Bohm scaling ($\delta \approx 1$), whereas in more quiescent downstream / cocoon regions (Zones B/C) we adopt a weaker energy dependence ($\delta \approx 1/2$ or $\delta \approx 1/3$). This zone dependence is motivated by recent discussions of microquasar-remnant transport and region-dependent diffusion [10]. The adopted coefficients are used consistently in the emissivity pipeline and clearly stated in the model setup.

2.2. Magnetic Field Prescription in Acceleration Regions

The magnetic field entering the acceleration and interaction rates is taken directly from the RMHD simulation output on a per-cell basis. To make the prescription explicit in acceleration zones, we use the following scheme:

• Primary: RMHD cell field. The local B-vector and magnitude are taken from the PLUTO RMHD solution in each cell.
• Equipartition/turbulent amplification control (acceleration cells). In cells tagged as acceleration sites (Zones A/B), we may enforce or verify that the magnetic energy density satisfies

  $$U_B={\displaystyle \frac{B^2}{8\pi }}={ϵ}_B{U}_{\mathrm{th}},$$

  with ${ϵ}_B$ a model parameter (typically $0.01$—1). If the RMHD cell field is below the chosen equipartition fraction in a flagged acceleration cell, we adopt a conservative floor $B\to max(B_{\mathrm{RMHD}},B_{\mathrm{eq}})$ with $B_{\mathrm{eq}}=\sqrt{8\pi {ϵ}_B{U}_{\mathrm{th}}}$. This makes the assumed field strength in acceleration regions explicit and reproducible.

At the jet base the initial toroidal field is $B_0={10}^4\mathrm{G}$, and the field evolves self-consistently in the RMHD simulation. The above prescription clarifies how B is interpreted/used in acceleration and emission calculations.

2.3. Nonthermal Proton Distribution

Nonthermal protons are described by a power-law energy distribution,

$$N_p(E)=K{n}_{p,0}{E}^{-\alpha },$$

with spectral index $\alpha \approx 2$ [3]. The normalization factor $K\ll 1$ specifies the fraction of hot protons relative to the thermal background density $n_{p,0}$. Proton energies extend up to $E_{max}={10}^6$ GeV.

Protons are assumed isotropic in the jet comoving frame when the scattering length is shorter than the radiative loss lengths [13]. Energy losses via mechanisms such as [3] synchrotron radiation, adiabatic expansion, and proton–proton collisions affect the steady-state proton spectrum.

Figure 1, Figure 2 and Figure 3 show proton density distributions and pion injection functions resulting from hadronic interactions.

2.4. Origin of Target Photon Fields

To remove ambiguity regarding the origin of photon fields used in $p\gamma$ interactions, we explicitly model the target photon field as the sum of physically motivated components:

• Companion star field (dominant in many HMXBs). We model the stellar radiation as a diluted blackbody with temperature $T_{\star }$ and radius $R_{\star }$. At a distance r from the star, the photon energy density scales as $u_{\star }(r)\propto {(R_{\star }/r)}^2$. The photon number density can be written as

  $$n_{\star }(ϵ,r)={\left({\displaystyle \frac{R_{\star }}{r}}\right)}^2{n}_{\mathrm{BB}}(ϵ,T_{\star }),$$

  where $n_{\mathrm{BB}}$ is the blackbody photon number density per unit energy.
• Accretion disk / corona component (when relevant). We include an additional thermal (or quasi-thermal) disk component with effective temperature $T_{\mathrm{disk}}$ and characteristic radius $R_{\mathrm{disk}}$, and optionally a coronal power-law tail if needed for higher-energy targets. In practice, this component can be switched on/off depending on the system class and parameter choice, and an accretions disk’s wind simplified construct is used instead.
• Scattered/wind photon field. A fraction of stellar/disk photons can be scattered in the wind environment, producing a more isotropized target field. We may treat this as a scaled component $n_{\mathrm{sc}}(ϵ)=f_{\mathrm{sc}}{n}_{\star }(ϵ,r)$ with $f_{\mathrm{sc}}\ll 1$.

For plasmoid–wind collisions, the stationary lab-frame photon distribution $n^{*}({ϵ}^{*})$ is Doppler-transformed into the jet comoving frame as described in Section 3.

3. Neutrino Emissivity

3.1. Neutrinos Produced Within the Jet

Neutrino production proceeds via pion production in proton–proton and proton–photon collisions, followed by pion decay. The proton distribution in each computational cell is Lorentz-transformed into the observer frame [16], and steady-state transport equations govern the pion and neutrino spectra [3,4,17,18].

The neutrino emissivity at neutrino energy E is given by

$$Q_{\pi \to \nu }(E)={\int }_E^{E_{max}}d{E}_{\pi }{t}_{\pi }^{-1}(E_{\pi })N_{\pi }(E_{\pi }){\displaystyle \frac{\Theta (1-r_{\pi }-x)}{E_{\pi }(1-r_{\pi })}},$$

where $N_{\pi }$ is the steady-state pion density, $t_{\pi }$ the pion decay timescale, $x=E/E_{\pi }$, $r_{\pi }={(m_{\mu }/m_{\pi })}^2$ and $\Theta$ the Heaviside function [4,19].

3.2. Neutrino Emissivity from Plasmoid–Wind Collision

To include neutrino emission enhancement from plasmoid collisions with ambient photon fields, the stationary photon distribution $n^{*}({ϵ}^{*})$ in the lab frame is Doppler-transformed into the jet comoving frame:

$$n(ϵ,\Omega )={\displaystyle \frac{n^{*}({ϵ}^{*})}{4\pi {\delta }_{*}^2}},$$

where ${\delta }_{*}$ is the Doppler factor calculated per cell using local velocity vector and line-of-sight angles, with the velocity reversed to represent incoming photons as seen in the jet frame [20]. The transformed photon distribution replaces the synchrotron photon field used in earlier approaches, and the target field is explicitly decomposed as in Section 2.4.

3.3. Target Photon Distribution and Doppler-Factor Cases Used in Figure 7

For the $p\gamma$ channel used in our SED calculations, we adopt a power-law target photon distribution in the lab/host frame and implement three prescriptions for the interaction Doppler factor $D_{\mathrm{ph}}$ (Cases a–c), as used in Figure 7. The target photon density entering the interaction kernel is written as

$$n_{\mathrm{ph}}(E_{\mathrm{ph}})=n_0{\left({\displaystyle \frac{E_{\mathrm{ph}}}{E_0}}\right)}^{-{\alpha }_{\mathrm{ph}}}{D}_{\mathrm{ph}}^{-2}.$$

(1)

Case (a): Constant Doppler factor. A single preset Doppler factor is assumed throughout the emitting plasmoid,

$$n_{\mathrm{ph}}^{(a)}=n_0{\left({\displaystyle \frac{E_{\mathrm{ph}}}{E_0}}\right)}^{-{\alpha }_{\mathrm{ph}}}{D}_{\mathrm{ph},0}^{-2}.$$

(2)

Case (b): Local Doppler factor (cell-based). The Doppler factor is computed locally for each cell using the cell velocity components and LOS angles,

$$n_{\mathrm{ph}}^{(b)}=n_0{\left({\displaystyle \frac{E_{\mathrm{ph}}}{E_0}}\right)}^{-{\alpha }_{\mathrm{ph}}}{\left[D(-u_{xx},-u_{yy},-u_{zz},{\varphi }_1,{\varphi }_2)\right]}^{-2}.$$

(3)

Case (c): Face-on collision approximation. The LOS angles are set to zero (head-on approximation),

$$n_{\mathrm{ph}}^{(c)}=n_0{\left({\displaystyle \frac{E_{\mathrm{ph}}}{E_0}}\right)}^{-{\alpha }_{\mathrm{ph}}}{\left[D(-u_{xx},-u_{yy},-u_{zz},0,0)\right]}^{-2}.$$

(4)

In all cases, the velocity vector is reversed in the Doppler-factor evaluation to represent photons entering the jet comoving frame (Appendix A).

3.4. Inclusion of Proton–Proton and Proton–Photon Contributions

Our emission modeling computes neutrino spectra including both proton–proton collisions (expected to be significant in matter-loaded/terminal regions) and proton–photon collisions with ambient photon fields.

Proton–proton pion injection follows

$$Q_{\pi }^{(pp)}(E_{\pi })=n_{p}c{\int }_{E_{\pi }}^{E_{p,max}}d{E}_p{N}_p(E_p){\sigma }_{pp}^{(\mathrm{inel})}(E_p)F_{\pi }^{(pp)}\left({\displaystyle \frac{E_{\pi }}{E_p}},E_p\right){\displaystyle \frac{1}{E_p}},$$

with ${\sigma }_{pp}^{(\mathrm{inel})}$ from [4] and pion distribution $F_{\pi }^{(pp)}$ as illustrated in Figure 2.

Proton–photon pion injection and neutrino production are computed analogously using the Doppler-transformed photon densities and the relevant cross sections and inelasticities [4]. In Section 6 we provide a direct comparison of the $pp$-only, $p\gamma$-only, and total neutrino outputs to quantify the role of each channel.

3.5. Velocity and Direction Filtering

To optimize computational costs, cells with velocity vectors aligned close to the observer’s line-of-sight and speeds exceeding $0.1c$ are prioritized for neutrino emission calculations. This filtering reduces the number of emitting cells, balancing fidelity and performance. A special “60 points case” restricts emission to a subset of points for benchmarking purposes.

3.6. Main Numerical Emission Expression (Implementation Form)

The explicit numerical form used in the emission code for the neutrino emissivity contribution arising from hadronic interactions (implemented in nemiss) can be written in a numerically convenient triple-integral form. Notation follows that of the main text: $E_{\nu }$ is the neutrino energy, $E_{\pi }$ the pion energy, $\xi$ the fraction $E_{\pi }/E_p$, $n_{\mathrm{dense}}^{\mathrm{fast}}$ the fast proton density function, $F_{\pi ,\mathrm{inj}}$ the pion injection kernel (single-collision pion distribution), ${\sigma }_{pp}^{\mathrm{inel}}$ the inelastic $pp$ cross-section, and $C_{Fbb\pi }$, $T_{\mathrm{depth}}^{pp}$ are implementation-specific normalization/depth factors used in the transport scheme. In the present implementation, the explicit form below corresponds to the $pp$ contribution; the $p\gamma$ contribution is computed analogously with the appropriate photon-target interaction kernel (Section and Appendix A).

$$\begin{array}{cc}\hfill q_{\nu }(E_{\nu };& u_{xx},u_{yy},u_{zz},{\varphi }_1,{\varphi }_2,b_{xx},b_{yy},b_{zz},n_{\mathrm{density}},p_{\mathrm{index}})=\hfill \\ \hfill & {\int }_{E_{\nu }}^{E_{max}}d{E}_{\pi }{\int }_{E_{\pi }/E_{p,max}}^{1}d\xi {\int }_{E_{p,min}}^{E_{p,max}}d{E}_p^{\prime }\Theta \left(1-r_{\pi }-{\displaystyle \frac{E_{\nu }}{E_{\pi }}}\right){\displaystyle \frac{t_{\pi }^{-1}(E_{\pi })}{E_{\pi }(1-r_{\pi })}}\hfill \\ \hfill & \times {\displaystyle \frac{c}{C_{Fbb\pi }(E_{\nu },b_{xx},b_{yy},b_{zz},n_{\mathrm{density}})}}{\displaystyle \frac{n_{\mathrm{dense}}^{\mathrm{fast}}\left(E_{\pi }/\xi ;p_{\mathrm{index}},u_{xx},u_{yy},u_{zz},{\varphi }_1,{\varphi }_2,n_{\mathrm{density}},b_{xx},b_{yy},b_{zz}\right)}{\xi }}\hfill \\ \hfill & \times F_{\pi ,\mathrm{inj}}\left(\xi ,E_{\pi }/\xi \right){\sigma }_{pp}^{\mathrm{inel}}\left(E_{\pi }/\xi \right)T_{\mathrm{depth}}^{pp}\left(E_{\nu },E_p^{\prime },E_{\pi },b_{xx},b_{yy},b_{zz},n_{\mathrm{density}}\right).\hfill \end{array}$$

(5)

This algebraic form is the implementation-level emissivity expression used in the code pipeline; it is algebraically equivalent to the standard cascade expression written in terms of $Q_{\pi }^{(pp)}$ and the pion decay kernel after the change of variables $\xi =E_{\pi }/E_p$ and inclusion of transport/depth normalizations. The detailed step-by-step derivation is provided in the Appendix.

4. Computer Programs Used

4.1. Rlos: Relativistic Line of Sight Imaging

The rlos code [21] performs special relativistic imaging by tracing rays through 4D RMHD simulation data, including relativistic beaming and time-delay effects.

4.2. PLUTO Hydrocode

PLUTO [22] is a shock-capturing, finite-volume RMHD code used here for the jet simulations on structured 3D meshes.

4.3. Nemiss

Nemiss [23,24] computes neutrino emissivities from hydrodynamic outputs, solving proton-to-neutrino cascades.

4.4. Additional Tools

Data visualization used Veusz. The codes are publicly available: PLUTO under GPL, nemiss and rlos under LGPL.

5. Model Setup

Our RMHD simulations model intermittent relativistic twin microquasar jets at $0.8c$. The jets propagate into ambient stellar and disk winds, with the companion star located outside the domain [25]. The initial magnetic field is toroidal with ${10}^4\mathrm{G}$ strength at the jet base.

Explicit model choices. The acceleration zones are defined as in Section 2.1. We use a zone-dependent diffusion parameterization $D(E)=D_0{(E/E_0)}^{\delta }$, adopting near-Bohm scaling ($\delta \simeq 1$) in plasmoid/shock regions and weaker scaling ($\delta \simeq 1/2$ or $1/3$) in downstream/cocoon regions (see Section 2.1). The magnetic field used in acceleration and interaction rates is taken from the RMHD cell field, with an optional equipartition floor in tagged acceleration cells (Section 2.2). Target photon fields, which may be explicitly decomposed, (Section 2.4) are transformed to the comoving frame where needed; for the $p\gamma$ SED cases, we additionally adopt the three Doppler prescriptions given in Section 3.3.

The computational mesh is $60\times 100\times 50$ Cartesian cells, each $2\times {10}^{10}$ cm in length. Simulation parameters are summarized in Table 1. Neutrino line-of-sight imaging uses local velocity and magnetic field data per cell, with relativistic effects modeled.

6. Results and Discussion

Intermittent twin jets propagate through ambient matter, generating dynamic equatorial structures as jet plasmoids interact with stellar and disk winds (Figure 4).

Synthetic neutrino images at $\sim 1.1$ GeV (Figure 5) reveal plasmoid-dominated emission enhanced by velocity filtering.

6.1. Relative Role of $pp$ and $p\gamma$ Channels

To quantify the significance of matter targets versus photon targets, we compare the neutrino outputs computed with $pp$-only, $p\gamma$-only, and combined ($pp+p\gamma$) emissivities. In matter-loaded and terminal regions, $pp$ can become competitive or dominant due to increased $n_p$, while in regions with intense external radiation fields the $p\gamma$ contribution can be enhanced.

Figure 6. Comparison of neutrino spectra computed with $pp$-only, $p\gamma$-only, and total ($pp+p\gamma$) emissivities at representative snapshots for the three collision scenarios shown in the top, middle, and bottom panels. All three panels are shown in unnormalized (relative) units and are intended solely to illustrate spectral shapes; the absolute normalization is introduced later in Figs. 7 and 10.).

Figure 6. Comparison of neutrino spectra computed with $pp$-only, $p\gamma$-only, and total ($pp+p\gamma$) emissivities at representative snapshots for the three collision scenarios shown in the top, middle, and bottom panels. All three panels are shown in unnormalized (relative) units and are intended solely to illustrate spectral shapes; the absolute normalization is introduced later in Figs. 7 and 10.).

6.2. Flux Normalization, Instrument Sensitivities, and Distance-Dependent Detectability

Figure 7 and Figure 8 present model-based neutrino emission results and their derived detectability. Figure 7 shows the predicted neutrino spectra for different emission prescriptions, while Figure 8 translates these spectra into an energy-dependent distance reach using instrument sensitivities. Figure 10 provides a complementary, instrument-only comparison in detector space, illustrating the relative normalized responses of current and future neutrino observatories across energy, independently of the specific model realization.

To assess detectability, the observable neutrino flux is obtained from the model luminosity via an assumed source distance d. Fluxes scale as

$${\Phi }_{\nu }(E;d)={\Phi }_{\nu }(E;d_0){\left({\displaystyle \frac{d_0}{d}}\right)}^2.$$

In figures that include observables, we overlay sensitivity curves for relevant instruments (e.g., IceCube, IceCube-Gen2, KM3NeT, Baikal-GVD) and we explicitly state the assumed distance used for flux conversion (here we use a canonical Galactic distance $d_0=8\mathrm{kpc}$ unless otherwise indicated).

A convenient distance-horizon metric at each energy is

$$d_{max}(E)=d_0\sqrt{{\displaystyle \frac{{\Phi }_{\nu }(E;d_0)}{{\Phi }_{\mathrm{sens}}(E)}}},$$

where ${\Phi }_{\mathrm{sens}}(E)$ is the instrument sensitivity at energy E. This provides a detectability-versus-distance statement independent of a single assumed distance.

Figure 9 shows the time evolution of total neutrino emissivity at 1.096 GeV.

Relativistic beaming enhances neutrino emission along local flow velocity directions, significantly affecting detectability for microquasars not perfectly aligned with Earth’s line of sight.

Figure 10. Normalized neutrino intensity as a function of energy for current and future neutrino telescopes (IceCube, IceCube-Gen2, KM3NeT, and Baikal-GVD). The red curve dashed shows model Case (b), normalized at the source to a total neutrino power $L_{\nu }={10}^{38}\mathrm{erg}{\mathrm{s}}^{-1}$ (corresponding to $\sim 1\%$ of the jet kinetic power), propagated to Earth assuming isotropic emission and mapped onto the normalized-intensity convention used for the experimental sensitivity curves. Shaded bands indicate the typical energy-dependent sensitivity envelopes of each instrument in normalized-intensity space, reflecting variations across analysis channels and effective area.

Figure 10. Normalized neutrino intensity as a function of energy for current and future neutrino telescopes (IceCube, IceCube-Gen2, KM3NeT, and Baikal-GVD). The red curve dashed shows model Case (b), normalized at the source to a total neutrino power $L_{\nu }={10}^{38}\mathrm{erg}{\mathrm{s}}^{-1}$ (corresponding to $\sim 1\%$ of the jet kinetic power), propagated to Earth assuming isotropic emission and mapped onto the normalized-intensity convention used for the experimental sensitivity curves. Shaded bands indicate the typical energy-dependent sensitivity envelopes of each instrument in normalized-intensity space, reflecting variations across analysis channels and effective area.

7. Conclusions

We present a time-resolved model for neutrino emission from relativistic microquasar jets combining RMHD simulations with explicit prescriptions for (i) acceleration zones and mechanisms, (ii) zone-dependent diffusion regimes, (iii) magnetic-field usage in acceleration regions, and (iv) physically motivated target photon fields. Neutrino production includes both proton–proton and proton–photon channels; their relative contributions are quantified, emphasizing that $pp$ can be significant in matter-loaded/terminal regions.

Synthetic neutrino images and spectra are assessed against the sensitivities of current and upcoming observatories, assuming a canonical Galactic distance and providing a distance-dependent detectability metric. Future work will explore higher resolution, refined radiation-field geometries, and more detailed system-specific parameter sets.

Appendix A.1. Proton–Proton Interaction and Pion Injection

$$Q_{\pi }^{(pp)}(E_{\pi })=n_{\mathrm{th}}c{\int }_{E_p^{\min }}^{E_p^{\max }}d{E}_p{N}_p(E_p){\sigma }_{pp}^{(\mathrm{inel})}(E_p)F_{\pi }^{(pp)}\left({\displaystyle \frac{E_{\pi }}{E_p}},E_p\right){\displaystyle \frac{1}{E_p}},$$

(A6)

where $n_{\mathrm{th}}$ is the density of thermal target protons, $N_p(E_p)$ is the differential number density of relativistic protons, ${\sigma }_{pp}^{(\mathrm{inel})}$ is the inelastic cross section, and $F_{\pi }^{(pp)}$ is the pion production spectrum per collision.

The inelastic cross section is approximated as [4]:

$${\sigma }_{pp}^{(\mathrm{inel})}(E_p)=\left(34.3+1.88L+0.25L^2\right){\left[1-{\left({\displaystyle \frac{E_{\mathrm{th}}}{E_p}}\right)}^4\right]}^2\times {10}^{-27}{\mathrm{cm}}^2,$$

(A7)

with $L=ln(E_p/1000\mathrm{GeV})$ and threshold energy $E_{\mathrm{th}}=1.22$ GeV.

Appendix A.2. Proton–Photon Interaction Rate

For the proton–photon channel, we use an isotropic photon field in the jet comoving frame. The interaction rate can be written as

$${\omega }_{p\gamma }(E_p)={\displaystyle \frac{c}{2{\gamma }_p^2}}{\int }_{{ϵ}_{\mathrm{th}}}^{\infty }d{ϵ}_r{\sigma }_{p\gamma }({ϵ}_r){\kappa }_{p\gamma }({ϵ}_r){ϵ}_r{\int }_{{ϵ}_r/(2{\gamma }_p)}^{\infty }dϵ{\displaystyle \frac{n_{\mathrm{ph}}(ϵ)}{{ϵ}^2}},$$

(A8)

where ${\gamma }_p=E_p/(m_p{c}^2)$, ${ϵ}_r$ is the photon energy in the proton rest frame, $n_{\mathrm{ph}}(ϵ)$ is the target photon density, ${\sigma }_{p\gamma }$ is the cross section, and ${\kappa }_{p\gamma }$ the inelasticity.

Appendix A.3. Doppler Transformations

Quantities computed in the comoving frame are transformed to the observer frame using the Doppler factor

$$\delta ={\displaystyle \frac{1}{\Gamma (1-\beta cos{\theta }_{\mathrm{obs}})}},$$

(A9)

where $\Gamma$ is the bulk Lorentz factor, $\beta =v/c$, and ${\theta }_{\mathrm{obs}}$ is the viewing angle.

The angle between the local velocity vector $\mathbf{v}$ and the line-of-sight unit vector $\mathbf{l}$ is:

$$cos\theta ={\displaystyle \frac{\mathbf{v}·\mathbf{l}}{|\mathbf{v}||\mathbf{l}|}}.$$

(A10)

The LOS vector components are:

$$l_x=cos{\varphi }_{2}cos{\varphi }_1,l_y=cos{\varphi }_{2}sin{\varphi }_1,l_z=sin{\varphi }_2.$$

(A11)

Photon energies transform as $E=\delta E^{\prime }$ and intensities as appropriate powers of $\delta$ depending on the specific observable discussed in the main text.

Appendix A.4. Line-of-Sight Geometry

The LOS unit vector $\mathbf{l}$ is parameterized as

$$l_x=cos{\varphi }_{2}cos{\varphi }_1,l_y=cos{\varphi }_{2}sin{\varphi }_1,l_z=sin{\varphi }_2.$$

The velocity vector is reversed in sign when calculating the Doppler factor for photons incoming into the jet comoving frame.

Appendix A.5. Some Model Setup Details

Appendix B.1. Proton–Proton Channel

We begin with the differential number density of relativistic protons, assumed to follow a power-law distribution in the jet comoving frame,

$$N_p(E_p)=K_p{E}_p^{-\alpha },$$

(A12)

where $K_p$ is the normalization constant and $\alpha$ the spectral index.

The pion injection rate from inelastic proton–proton collisions is

$$Q_{\pi }^{(pp)}(E_{\pi })=n_{\mathrm{th}}c{\int }_{E_{\pi }}^{E_{p,max}}d{E}_p{N}_p(E_p){\sigma }_{pp}^{\mathrm{inel}}(E_p)F_{\pi }^{(pp)}\left({\displaystyle \frac{E_{\pi }}{E_p}},E_p\right){\displaystyle \frac{1}{E_p}},$$

(A13)

where $n_{\mathrm{th}}$ is the density of thermal target protons and $F_{\pi }^{(pp)}$ is the pion production spectrum per collision.

Assuming steady state, the pion distribution follows from injection–decay balance,

$$N_{\pi }(E_{\pi })={\displaystyle \frac{Q_{\pi }^{(pp)}(E_{\pi })}{t_{\pi }^{-1}(E_{\pi })}},$$

(A14)

where $t_{\pi }$ is the pion decay timescale.

The neutrino emissivity from charged pion decay is

$$Q_{\nu }^{(pp)}(E_{\nu })={\int }_{E_{\nu }}^{E_{\pi ,max}}d{E}_{\pi }{\displaystyle \frac{t_{\pi }^{-1}(E_{\pi })N_{\pi }(E_{\pi })}{E_{\pi }(1-r_{\pi })}}\Theta \left(1-r_{\pi }-{\displaystyle \frac{E_{\nu }}{E_{\pi }}}\right),$$

(A15)

with $r_{\pi }={(m_{\mu }/m_{\pi })}^2$. Substituting $N_{\pi }(E_{\pi })$ yields the emissivity integral implemented numerically.

Appendix B.2. Proton–Photon Channel

For the proton–photon channel, the interaction rate in an isotropic photon field is

$${\omega }_{p\gamma }(E_p)={\displaystyle \frac{c}{2{\gamma }_p^2}}{\int }_{{ϵ}_{\mathrm{th}}}^{\infty }d{ϵ}_r{\sigma }_{p\gamma }({ϵ}_r){\kappa }_{p\gamma }({ϵ}_r){ϵ}_r{\int }_{{ϵ}_r/(2{\gamma }_p)}^{\infty }dϵ{\displaystyle \frac{n_{\mathrm{ph}}(ϵ)}{{ϵ}^2}},$$

(A16)

where ${\gamma }_p=E_p/(m_p{c}^2)$, ${ϵ}_r$ is the photon energy in the proton rest frame, and $n_{\mathrm{ph}}(ϵ)$ is the target photon density.

The pion injection rate is then

$$Q_{\pi }^{(p\gamma )}(E_{\pi })={\int }_{E_{\pi }}^{E_{p,max}}d{E}_p{N}_p(E_p){\omega }_{p\gamma }(E_p)\delta \left(E_{\pi }-{\kappa }_{\pi }{E}_p\right),$$

(A17)

where ${\kappa }_{\pi }$ is the average inelasticity mapping proton to pion energy.

Following the same steady-state treatment as in the $pp$ case, the neutrino emissivity is

$$Q_{\nu }^{(p\gamma )}(E_{\nu })={\int }_{E_{\nu }}^{E_{\pi ,max}}d{E}_{\pi }{\displaystyle \frac{t_{\pi }^{-1}(E_{\pi })N_{\pi }(E_{\pi })}{E_{\pi }(1-r_{\pi })}}\Theta \left(1-r_{\pi }-{\displaystyle \frac{E_{\nu }}{E_{\pi }}}\right).$$

(A18)

Appendix B.3. Consistency with Numerical Implementation

Both channels reduce to the same formal structure for neutrino production, differing only in the pion injection term. The expressions above correspond directly to the numerical emissivity kernels implemented in the simulation pipeline and ensure energy conservation and consistency with standard hadronic cascade theory.

Appendix D.1. Source-Level Normalization

The neutrino spectra presented in Fig. 7 are interpreted as differential neutrino power distributions per logarithmic energy interval,

$${\mathcal{P}}_{\nu }(E)\equiv E^2{\displaystyle \frac{d{\dot{N}}_{\nu }}{dE}},$$

(A19)

with units of $\mathrm{erg}{\mathrm{s}}^{-1}$. We assume that a fraction ${\eta }_{\nu }$ of the jet kinetic power $L_{\mathrm{jet}}$ is converted into neutrinos, so that the total neutrino luminosity is

$$L_{\nu }={\eta }_{\nu }{L}_{\mathrm{jet}}.$$

(A20)

The model spectrum is normalized by enforcing

$${\int }_{E_{min}}^{E_{max}}{\mathcal{P}}_{\nu }(E){\displaystyle \frac{dE}{E}}=L_{\nu },$$

(A21)

i.e., the integral over $dlnE$ equals the assumed neutrino power.

Appendix D.2. Flux at Earth and Comparison Units

For a source at distance d, the corresponding spectral energy flux is

$$E^2{\displaystyle \frac{d{\Phi }_{\nu }}{dE}}={\displaystyle \frac{{\mathcal{P}}_{\nu }(E)}{4\pi d^2}},$$

(A22)

which preserves the normalization in physical units.

To compare with the normalized-instrument presentation of Fig. 10, we map the model curve into the same plotting scale by applying a single multiplicative factor C (defined in the main text / figure caption) so that

$$I_{\nu }^{(\mathrm{norm})}(E)=C\left(E^2{\displaystyle \frac{d{\Phi }_{\nu }}{dE}}\right),$$

(A23)

and the resulting curve is overlaid against the instrument sensitivity bands without altering its intrinsic spectral shape.

Appendix D.3. Case Used

Unless otherwise stated, the normalization shown in Fig. 10 corresponds to model case (b) as discussed in the main text.