© G. J. Escobar et al. 2022

1. Introduction

The origin of Galactic cosmic rays (CRs) is still a matter of debate. At present, supernova remnants (SNRs) are considered the most plausible sites where these particles accelerate up to high energies that then propagate diffusively through the ambient medium. Even though this paradigm explains the CR power observed from Earth, it fails to simultaneously reproduce the main characteristics observed in the CR spectrum, such as the energy at which the transition to extragalactic CRs is found, the slope of the spectrum, or the observed anisotropy and chemical composition (Ginzburg & Syrovatskii 1963; see also Blasi 2013 for a recent review). This lack of agreement may reflect the existence of alternative sources of CRs, but also a poor knowledge of either the dynamics of CRs or the properties of the ambient medium through which these particles propagate.

Various attempts have been made to explain the observed CR spectrum based on contributions from different sources of particle acceleration. Proposals include that CRs originate in pulsars (e.g., Bednarek 2002; Bednarek & Bartosik 2004), microquasars (MQs; Hillas 1984; Heinz & Sunyaev 2002; Bednarek 2005), winds of massive stars or stellar clusters, and enhanced star-forming regions (e.g., Morlino et al. 2021; Peretti et al. 2022). Romero & Vila (2009) also suggested that MQs might be responsible for the positron component of the CR spectrum. These scenarios are based on the presence of outflows developing strong magnetized shocks in which particles accelerate to relativistic energies. The same magnetic fields are, however, also efficient in trapping and cooling particles, therefore reducing the power and maximum particle energies achieved.

Recently, we developed a different scenario (Escobar et al. 2021) in which the escape of CRs from magnetically confined jets of MQs is mediated by relativistic neutrons. The presence of charged particles accelerated to relativistic energies in MQ jets follows from observations of the non-thermal radio and gamma-ray emission of these systems. Observations show that hadronic content is present at least in some MQs (Migliari et al. 2002; Díaz Trigo et al. 2013), which is also supported by theoretical models of jet launching mechanisms (Blandford & Payne 1982). The production of relativistic neutrons is an inescapable consequence of these facts, given that a fraction of the interactions between relativistic and thermal protons produce neutrons. These neutrons carry most of the energy of the progenitor projectiles and escape freely from magnetic fields that confine the jet. Indeed, relativistic neutron production up to ∼EeV energies has been explored in the jets of active galactic nuclei (AGNs; e.g., Sikora et al. 1989; Atoyan & Dermer 2003), although the large spatial scales of these systems prevent that neutrons escape before decaying. Our model predicts that MQ jets with large kinetic luminosities and low velocities produce neutrons up to PeV energies that decay into protons and electrons outside the system, therefore directly injecting CRs in the surrounding medium (Escobar et al. 2021). These CRs display steep spectra (i.e. spectral index ≳2.5), and may carry up to ∼10^{48 − 49} erg throughout the whole MQ life, only ∼1 − 2 orders of magnitude less than a SNR.

The scenario devised by Escobar et al. (2021) is based on the MQ jet model of Romero & Vila (2008). Its main ingredient is inelastic p–p scattering, which produces relativistic neutrons from the abundant relativistic protons accelerated in jet internal shocks. This process competes mainly with proton adiabatic cooling (except at sites near the jet base and very high energies at which proton synchrotron dominates cooling) that is assumed to operate in a conic jet. Some features of the model are adopted from the observed properties of one of the best studied MQs, namely Cygnus X-1 (Pepe et al. 2015). Recently, Kantzas et al. (2021) developed a different model for the same MQ, which includes some ingredients not present in that of Romero & Vila (2008). The first one presents a parabolic geometry for the jet, which directly affects the adiabatic cooling regulating losses and should therefore affect neutron production. Also, convection within the jet is treated in the model of Kantzas et al. (2021). This may transport relativistic protons to lower density regions, decreasing the p–p scattering rate. Finally, second-order processes, such as photomeson production in the synchrotron radiation field of the jet, have been shown to contribute to γ-ray emission by these authors and may therefore contribute also to the neutron production. These processes are absent from the work of Escobar et al. (2021), where only estimates of the contribution from external photon fields are implemented. In this work, we improve upon the model devised by these authors in order to explore the effects of the jet geometry, convection, and p–γ interactions on CR production mediated by relativistic neutrons.

This paper is organised as follows. In Sect. 2, we describe the jet model and the neutron production mechanisms. In Sect. 3, we explore the neutron and cosmic-ray luminosities and spectra for different scenarios in our model. We discuss the results and present our conclusions in Sect. 4.

2. Jet model

We used a MQ jet model based on the one built and described in Escobar et al. (2021). We refer to this work for details of the original model, as here we focus on the new developments. We considered a lepto-hadronic jet of luminosity L_jet and bulk velocity v_jet corresponding to a bulk Lorentz factor Γ. The jet is launched at a distance z₀ from the compact object (see Fig. 1). The magnetic field, B, is computed assuming equipartition between kinetic and magnetic energies at the jet base (z₀), and then adopting a power-law function of the distance along the jet, B ∝ z^−m, with m serving as the magnetic index. A fraction of the jet power is deposited in relativistic particles once the flow develops an acceleration mechanism. The distance at which this acceleration region begins is denoted by z_a, and its top by z_t, while the emission region extends up to z_max ≥ z_t. The shape of the emission region is parameterised by

$$\begin{array}{c}\hfill r(z)=r_0{\left(\frac{z}{z_0}\right)}^{\alpha },\end{array}$$(1)

Fig. 1. Schematic picture of the jet (not to scale), showing a jet launched at a distance z0 from the compact object. The acceleration region, shaded in red, is defined between za and zt, while the emission region spans over the range za ≤ z ≤ zmax.

where r₀ is the jet radius at its base and 0 < α ≤ 1 describes the jet geometry. The jet collimation increases with decreasing α; the values α = 1/2 and α = 1 represent parabolic and conic shapes, respectively, and the jet becomes almost cylindrical as α → 0. This prescription allows us to vary the strength of adiabatic cooling, which competes with the p–p interactions that produce neutrons.

We assume that the relativistic population takes 10% of the jet luminosity (q_rel = 0.1). This fraction is distributed between protons and electrons according to a hadron-to-lepton power-ratio a. The discussion about the acceleration mechanism is beyond the scope of this work; here, we only mention that particles accelerated up to relativistic energies in astrophysical sources are likely to develop a power-law distribution once a steady state is achieved (e.g., Schlickeiser 2002). Therefore, assuming that the acceleration time scale is much shorter than those of other processes, we parameterise the injection as a power law in the proton energy, E_p, with spectral index p and an exponential cutoff at energy, $E_{\mathrm{p}}^{\max }(z)$:

$$\begin{array}{c}\hfill Q_{\mathrm{p}}(E,z)=Q_0{E}_{\mathrm{p}}^{-p}exp(-\frac{E_{\mathrm{p}}}{E_{\mathrm{p}}^{\max }(z)})f(z).\end{array}$$(2)

The maximum energy $E_{\mathrm{p}}^{\max }$ is achieved when the total loss rate equals the acceleration rate of the particles (see discussion below), and Q₀ is obtained by normalising the injection function to the total power of the relativistic proton population. The function f(z) is chosen to represent the result of the convective motion that produces particle injection beyond the acceleration region. Therefore, we chose f(z) = 1 for z_a ≤ z ≤ z_t and a power-law f(z) = (z/z_t)^−δ for z > z_t. The index δ > 0 is left free to vary the convective strength (lower values representing stronger convection; see discussion in Sect. 3).

The steady-state solution for the densities of the relativistic hadron populations is given by the coupled transport equations:

$$\begin{array}{cc}\hfill \frac{\partial }{\partial E}[b_{\mathrm{p}}{N}_{\mathrm{p}}]+t_{\mathrm{esc}}^{-1}{N}_{\mathrm{p}}& =Q_{\mathrm{p}}-{\mathrm{\Lambda }}_{\mathrm{p}},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill t_{\mathrm{esc},\mathrm{n}}^{-1}{N}_{\mathrm{n}}& =Q_{\mathrm{n}}.\hfill \end{array}$$(4)

Here, N_p(n) is the proton (neutron) steady-state particle density, t_esc and t_esc, n are the proton and neutron escape rates, respectively, and b_p accounts for the total energy loss rate for protons. The terms Q_n and Λ_p represent the source rate of neutrons and the related sink rate of protons, respectively. All quantities depend on the energy of the corresponding particle population and the distance, z. The boundary condition is $N_{\text{p}}(E_{\text{p}}^{\text{max}}(z),z)=0$. The effect of any term depending on the derivatives with respect to z is included effectively by the function f(z) and therefore not included explicitly in Eq. (4).

Proton-proton inelastic collisions are a source mechanism for relativistic neutrons through the branching channels in which the relativistic proton is destroyed, yielding a neutron that carries approximately half of the proton energy plus several secondary particles. The total cooling rate is:

$$\begin{array}{c}\hfill t_{\mathrm{pp}}^{-1}={\sigma }_{\mathrm{pp}}{n}_{\mathrm{p}}{v}_{\mathrm{rel}}{K}_{\mathrm{pp}},\end{array}$$(5)

where σ_pp and K_pp are the cross section and inelasticity of the process, v_rel is the relative velocity of colliding particles, which we take as c, and then n_p = n_b + n_w is the target proton density, which accounts for the contribution of protons from the jet bulk material (n_b) and the wind of the companion star (n_w). These densities are given by:

$$\begin{array}{cc}& n_{\mathrm{b}}=\frac{L_{\mathrm{jet}}(1-q_{\mathrm{jet}})}{m_{\mathrm{p}}{c}^2{v}_{\mathrm{jet}}\pi r{(z)}^2},\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill & n_{\mathrm{w}}=\frac{\mathrm{\Gamma }\dot{M}}{4\pi (d^2+z^2)v_{\mathrm{w}}{m}_{\mathrm{H}}},\hfill \end{array}$$(7)

in the jet reference frame. In the latter equation, Ṁ is the mass-loss rate of the companion, v_w its wind velocity, d its distance to the compact object, and m_H is the mass of the hydrogen atom. We adopted the standard velocity profile for a line-driven wind given by Lamers & Cassinelli (1999). A discussion about the impact of the wind density on the model and its results can be found in Escobar et al. (2021). Both contributions are treated separately to assess their individual significances. We adopted the parameterisation of the cross section given by Kafexhiu et al. (2014). Regarding the neutron channel, we adopted a conservative branching ratio of 0.16 (see discussion in Sect. 2 of Escobar et al. 2021). Thus, the proton sink rate due to this process is given by ${\mathrm{\Lambda }}_{\mathrm{pp}}(E_{\mathrm{p}})=0.16t_{\mathrm{pp}}^{-1}(E_{\mathrm{p}})N_{\mathrm{p}}(E_{\mathrm{p}})$ and the corresponding neutron source rate is $Q_{\mathrm{n}}^{(\mathrm{pp})}(E_{\mathrm{n}})$ = Λ_pp(E_p), with E_n = 0.5E_p.

Some channels of photohadronic interactions also produce relativistic neutrons. We considered the inelastic collisions between protons and radiation fields, both internal and external to the jet. To compute the cooling rate, inelasticities, and inclusive cross sections, we adopted the prescription given in Atoyan & Dermer (2003). For the single-pion channel (photon energies ε between 200 MeV and 500 MeV, in the reference frame of the proton), the cross-section can be approximated as σ₁ ≈ 340 μbarn and the inelasticity as K₁ = 0.2. In the multi-pion channel (ε ≥ 500 MeV), these values are σ₂ ≈ 120 μbarn and K₂ ≈ 0.6. For both channels, the probability of conversion of a proton to a neutron is ≈0.5. The cooling rate is given by:

$$\begin{array}{c}\hfill t_{\mathrm{p}\gamma }^{-1}(E_{\mathrm{p}})={\displaystyle {\int }_{\frac{{\epsilon }_{\mathrm{th}}}{2{\gamma }_{\mathrm{p}}}}^{\infty }\mathrm{d}\epsilon \frac{n_{\mathrm{ph}}c}{2{\epsilon }^2{\gamma }_{\mathrm{p}}^2}{\int }_{{\epsilon }_{\mathrm{th}}}^{2\epsilon {\gamma }_{\mathrm{p}}}{\epsilon }^{\prime }{\sigma }_{\mathrm{p}\gamma }({\epsilon }^{\prime })K_{\mathrm{p}\gamma }({\epsilon }^{\prime })\mathrm{d}{\epsilon }^{\prime },}\end{array}$$(8)

where γ_p is the proton Lorentz factor, ε the target photon energy, ε_th the energy threshold, n_ph the differential number density of target photons, and σ_pγ and K_pγ are the cross section and inelasticity functions, respectively. Similarly, the collision rate is given by:

$$\begin{array}{c}\hfill {\nu }_{\mathrm{p}\gamma }(E_{\mathrm{p}})={\displaystyle {\int }_{\frac{{\epsilon }_{\mathrm{th}}}{2{\gamma }_{\mathrm{p}}}}^{\infty }\mathrm{d}\epsilon \frac{n_{\mathrm{ph}}c}{2{\epsilon }^2{\gamma }_{\mathrm{p}}^2}{\int }_{{\epsilon }_{\mathrm{th}}}^{2\epsilon {\gamma }_{\mathrm{p}}}{\epsilon }^{\prime }{\sigma }_{\mathrm{p}\gamma }({\epsilon }^{\prime })\mathrm{d}{\epsilon }^{\prime }.}\end{array}$$(9)

The proton and neutron energies are related by $E_{\mathrm{n}}=(1-{\overline{K}}_{\mathrm{p}\gamma })E_{\mathrm{p}}$, where ${\overline{K}}_{\mathrm{p}\gamma }={\nu }_{\mathrm{p}\gamma }^{-1}{t}_{\mathrm{p}\gamma }^{-1}$ is the mean inelasticity. The proton sink rate of this interaction is given by ${\mathrm{\Lambda }}_{\mathrm{p}\gamma }=0.5t_{\mathrm{p}\gamma }^{-1}{N}_{\mathrm{p}}$, while the neutron source rate is $Q_{\mathrm{n}}^{(\mathrm{p}\gamma )}(E_{\mathrm{n}})={\mathrm{\Lambda }}_{\mathrm{p}\gamma }(E_{\mathrm{p}})$. Finally, the total proton sink rate and total neutron source rate are given by Λ_p = Λ_pp + Λ_pγ and $Q_{\mathrm{n}}=Q_{\mathrm{n}}^{(\mathrm{pp})}+Q_{\mathrm{n}}^{(\mathrm{p}\gamma )}$, respectively.

The relativistic populations also suffer non-radiative losses, such as escape, decays, and adiabatic work. We estimated the escape rate of protons as $t_{\mathrm{esc}}^{-1}=v_{\mathrm{jet}}{(z_{\max }-z)}^{-1}$, representing the loss through the head of the jet. The fraction of particles that escape sideways is of the order of r_g/r(z)≪1, with r_g the gyroradius, so we neglect it. The previous condition is not fulfilled by particles with energies close to the maximum energy; the power carried by these particles, however, is negligible with respect to the total power of the population. We modelled the escape rate of neutrons as $t_{\mathrm{esc},\mathrm{n}}^{-1}=c{(z_{\max }-z)}^{-1}$. As disussed in Escobar et al. (2021), this is a conservative estimate, even though our main results do not depend on the exact value of such quantity. In addition, we neglected any neutron decay within the jet due to the small size of the latter compared to the mean decay length at any relevant energy. On the other hand, the adiabatic loss-rate for protons is expressed as:

$$\begin{array}{c}\hfill t_{\mathrm{ad}}^{-1}=\frac{2\alpha }{3z}{v}_{\mathrm{jet}},\end{array}$$(10)

given that these particles are relativistic.

The acceleration rate can be written as

$$\begin{array}{c}\hfill t_{\mathrm{acc}}^{-1}=\frac{\eta ceB}{E},\end{array}$$(11)

where η is an efficiency parameter and e is the elementary charge. Protons also cool through radiative losses, including the already mentioned inelastic collisions with matter and radiation, and synchrotron emission. The formulae for the synchrotron losses can be found in Blumenthal & Gould (1970). Overall, neutrons suffer neither adiabatic nor radiative losses.

3. Neutron and cosmic-ray production

In this section, we investigate the production of relativistic neutrons for several scenarios comprising different features of the emission region, the convection of particles, the collimation of the jet, and the presence of different photon fields. Table 1 shows the parameters for the fiducial case, which is nearly the same used in Escobar et al. (2021), except for the minimum energy and acceleration efficiency, whose values were adopted to ensure that in all cases acceleration overtakes losses. This scenario has been adapted from Pepe et al. (2015) to represent the jet of Cygnus X-1. Regarding the features under investigation, this is a conic jet with a compact acceleration region near its base, which is also the emission region (i.e. convection is negligible). No photon fields are taken into account in this case; we include the contribution of photohadronic interactions separately in Sect. 3.3.

3.1. Acceleration region and convection

We first explored several models with different sizes for the acceleration region, starting from a compact one with z_t = 1.9 × 10¹² cm, and then extending it up to z_t = 10¹⁵ cm. The last value is the extension of the jet of Cygnus X-1 reported by Stirling et al. (2001), showing that particle acceleration might develop at much higher distances than those reported by Pepe et al. (2015; cf. Kantzas et al. 2021). In Fig. 2, we show the proton loss and acceleration rates at different sites of the jet for the model with the most extended acceleration region. It is important to stress that except for escape, energy losses at a given position z in the jet do not vary with the size of the emission region, therefore, the loss curves of Fig. 2 apply to any of the aforementioned models.

Fig. 2. Proton cooling rates for the case zt = zmax = 1015 cm, at the base (left panel), logarithmic midpoint (middle panel), and the end of the acceleration region (right panel).

We observe that inelastic p–p losses, responsible for the production of relativistic neutrons, are larger close to the base of the jet and decrease upwards. This is due to the decreasing density of bulk protons as the jet expands. The most important process that competes against p–p is adiabatic energy loss, except for the highest energies at the jet base, where synchrotron radiation dominates, and for all energies at the jet top, where the escape losses slightly exceed adiabatic ones. Adiabatic loss and p–p scattering are of the same order of magnitude at the jet base, but upwards p–p losses decrease faster than adiabatic work, becoming orders of magnitude smaller. These features suggest that conditions for neutron production are most favourable at the base.

Figure 3 supports the interpretation expressed above, showing that the fraction of power per unit volume injected relativistic protons, q_p, that is transferred to neutron emissivity, q_n, decreases with z, where:

$$\begin{array}{c}\hfill q_{\mathrm{p}(\mathrm{n})}(z)={\displaystyle \int }{E}_{\mathrm{p}(\mathrm{n})}{Q}_{\mathrm{p}(\mathrm{n})}(E_{\mathrm{p}(\mathrm{n})},z)\mathrm{d}{E}_{\mathrm{p}(\mathrm{n})}.\end{array}$$(12)

Fig. 3. Ratio of neutron emissivity to proton volumetric injection rate (top) and neutron luminosity distribution normalised with total luminosity Ln (bottom).

Despite the fact that the jet base is more efficient in transferring energy to the neutron channel, the total neutron luminosity, L_n, is still dominated by the jet top due to its far larger volume. At the same total power injected into relativistic particles, these facts render extended acceleration regions less efficient than compact ones for neutron production. This is because in compact acceleration regions, most of the power is released near the jet base, where neutrons can be efficiently produced; whereas extended ones inject most of their energy in large regions far from the base, where neutron production is not efficient. In Fig. 4, we show that indeed the fraction of the total luminosity of relativistic particles transferred into neutrons decreases with increasing acceleration region size. The linear trend with log z_t arises from the fact that for a conic jet adiabatic losses dominate the proton cooling (no dependence with z); for neutrons, the dominant loss process is escape (varies according to ∼z⁻¹). Therefore, jets with compact acceleration regions near their base are the best candidates for injecting neutrons into their surroundings.

Fig. 4. Fraction of the luminosity of relativistic protons transferred to neutrons, as a function of the size of the acceleration region, defined by the position of its top zt. Here, Lp stands for the total luminosity in the proton relativistic population.

It is important to recall that this behaviour is driven by the competition between p–p scattering and adiabatic losses. The latter are directly related to the jet degree of collimation (see Eq. (10)) and therefore we expect that collimation plays a major role in neutron production. We defer this discussion to Sect. 3.2

The top panel of Fig. 5 shows the maximum energy $E_{\text{p}}^{\text{max}}$ achieved by protons as a function of the position, z, for three models with different sizes of the acceleration region. The behaviour of this function reflects the existence of three regimes in which different loss processes dominate: near the jet base (z ≲ 10⁹ cm) $E_{\text{p}}^{\text{max}}$ increases with z because synchrotron losses dominate at the higher proton energies while the magnetic field drops with z. At further distances, the energy loss is ruled mainly by adiabatic work, causing $E_{\text{p}}^{\text{max}}$ to decrease with z with a constant slope (z ≈ 10⁹ cm to z ≈ 10¹⁴ cm). Finally, near the top of the acceleration region, the escape rate overcomes the adiabatic losses and limits $E_{\text{p}}^{\text{max}}$, steepening its decrease (z ≳ 10¹⁴ cm). The distance z at which the slope breaks increases with z_t, because escape depends on the size of the emission region (and z_max = z_t in this discussion). More compact regions have higher escape rates, resulting in lower maximum proton energies. Based on this behaviour, we expect the neutron spectrum to be very steep and dominated by low-energy neutrons produced far from the jet base in the general case, reaching high energies only for compact acceleration regions.

Fig. 5. Maximum energy of protons as a function of distance z along the jet for models comprising three different values of the acceleration region top zt. The top (bottom) panel corresponds to a magnetic index m = 1.9 (1.5).

The top panel of Fig. 6 shows the integrated neutron spectra of the jet for different values of the top of acceleration region, z_t (=z_max in this scenario). In this case, we adopt E_min = 5 m_pc². Neutron spectra have common features independently of the extension of the emission region: a broken power-law shape with similar spectral indices among all the different scenarios. An index of 3.4, indicating a very steep behaviour at high energies as predicted by our previous discussion, as well as a shallower one (1.3) at low energies, describe extremely well all curves. These features are remarkably independent of the characteristics of the scenario. The break in the spectrum is related to the minimum value of $E_{\text{p}}^{\text{max}}$ in each model, which is also $E_{\text{p}}^{\text{max}}(z_t)$ and therefore indicates the shrinking of the effective emission region. Neutrons of lower energies are produced in the whole emission region, whereas those of higher energies cannot be produced in the external parts having lower $E_{\text{p}}^{\text{max}}$. Below the break, the power emitted in neutrons decreases with energy only due to the corresponding change in acceleration and loss mechanisms. The low-energy index is indeed directly related to the power-law slope of the injected proton energy spectrum due to the weak dependence of losses on the proton energy (1.3 ≈ 1.4 = p − 1, the decrement in one arising from a factor E_p between the number and the energy of protons). Above the break, the change in the effective emission volume produces a faster decrease. The high-energy index is therefore determined by the interplay between the change with z of the maximum proton energy and that of the cross section of the jet.

Fig. 6. Integrated neutron power spectra in the jet frame, for different values of the top of the acceleration (and in this case also emission) region, zt = zmax. Blue stars indicate the value of Pn evaluated at the maximum neutron energy at zt = zmax, i.e. the minimum value of $E_{\mathrm{n}}^{\max }(=0.5E_{\mathrm{p}}^{\max })$ along the jet for each model (cf. Fig. 5). In all the cases Emin = 5mc2. The top (bottom) panel corresponds to a magnetic index m = 1.9 (1.5). The remaining parameters are those of Table 1.

The above reasoning implies that the value of the high-energy index of the neutron spectrum is ultimately governed by the structure of the jet magnetic field, that determines the variation of $E_{\text{p}}^{\text{max}}$. The bottom panel of Fig. 5 shows the maximum energy of protons along the jet for a different value of the magnetic index, m = 1.5. A strong difference is seen with respect to the case of m = 1.9. A weaker decrease of B with z, expressed by a lower value of m, translates into a correspondingly weaker decrease of $E_{\text{p}}^{\text{max}}$. This enlarges the energy range at which the full jet volume is emitted, moving the spectral break to higher neutron energies (bottom panel of Fig. 6). Below the break, the spectral index remains unchanged for different values of the magnetic index, as the former depends mainly on the injection index p. Since the total proton luminosity is fixed, the enlargement of the full-volume spectral range leaves less power available above the break, steepening the high-energy part of the spectrum.

We also studied the effect of different convection regimes on neutron production. For this purpose, we used models with a compact acceleration region with z_t = 1.9 × 10¹² cm and an emission region extending to z_max = 10¹⁵ cm. Convection transports relativistic particles from the former to the latter. As convection is stronger, more power in relativistic protons is transported to regions far from the base, where the density of bulk protons is lower and therefore p–p collisions are less frequent. Thus, we expect stronger convective jets to be less efficient producing neutrons. In Fig. 7, we show the total power of produced neutrons, normalised to the injected proton power, for different convection regimes. In fact, the plot shows that a weakly convective behaviour is more efficient in producing neutrons that a strong one, the power released in neutrons being about three orders of magnitude higher in the former case.

Fig. 7. Fraction of the luminosity of relativistic protons transferred to neutrons, as a function of the convective strength of the jet, defined by the parameter δ.

In Fig. 8, we show the neutron spectra comprising the different convection scenarios. As in the previous case, we use E_min = 5 m_pc² (and with the remaining properties of the system as in the fiducial case). Comparing with the fiducial scenario in Fig. 6, we see that at values of δ ≈ 4 convection can be considered negligible. All neutron spectra show the same common features found in Fig. 6: a broken power-law shape with similar spectral indices among the different models. The spectral indices are independent of the convection strength, and their values agree with those of the previous scenario with no convection. The power law breaks at higher energies in the case of a weaker convection regime, which is consistent with our previous finding, because a weaker convection implies a smaller effective emission region.

Fig. 8. Integrated neutron power spectra in the jet frame, for scenarios with different convective strength, with Emin = 5mc2. The remaining parameters are the same as those in Table 1.

3.2. Collimation

We go on to analyse the effect of jet collimation in the production of neutrons. We mention above that the neutron power results mainly from the competition between p–p collisions and adiabatic losses. The latter depend strongly on the jet collimation and if they are reduced, it is expected that proton-proton collisions dominate over a wide range of distances near the base of the acceleration region, increasing the power transferred to neutrons. Escape provides another, but less relevant, competing process. We note from Fig. 2 that the rates of these three processes are roughly constant along the whole range of proton energies. Thus, such losses are well represented by an average rate over energies, which we denote by ${\overline{t}}^{-1}$. In order to measure the relevance of proton-proton collisions with respect to other losses for different degrees of jet collimation, we compute the following ratio:

$$\begin{array}{c}\hfill f=\frac{{\overline{t}}_{\mathrm{pp}}^{-1}}{{\overline{t}}_{\mathrm{ad}}^{-1}+{\overline{t}}_{\mathrm{esc}}^{-1}}.\end{array}$$(13)

Synchrotron losses are not taken into account since they are not relevant for this analysis, because they only contribute very near the jet base at the highest energies. Figure 9 shows f as a function of the distance along the jet, z, for the scenario where the jet varies from low (α = 1) to high collimation (α → 0) through the shape parameter α. Also, different bulk Lorentz factors are explored, as they also determine the level of adiabatic losses. The range is the same as in Escobar et al. (2021), Γ_min = 1.034 (e.g., that of the jet in the microquasar SS 433; Chaty 2007) and Γ_max = 5. The remaining parameters are those of the fiducial case. It can be seen that independently of the value of α, higher p–p loss rates are obtained if Γ is lower; this is in accordance with results of the previous work (Escobar et al. 2021). The shape of the jet has also a major impact on the energetic loss balance. In all cases, the p–p loss rate is more relevant at regions closer to the jet base, where matter target fields are denser. Well-collimated jets reduce adiabatic losses, redirecting energy to the p–p channel. In the case α = 0.07 (an almost cylindrical jet), proton-proton collisions dominate over the rest of the losses for almost all values of Γ and over the whole region of acceleration. This suggests that highly collimated jets are efficient neutron production sites, which is confirmed by Fig. 10. We observe that in the case of a parabolic jet, the total power in the neutron population rises by about four orders of magnitude with respect to the conic case.

Fig. 9. Ratio of proton–proton to adiabatic-plus-escape average loss rates, f, for models with different collimation parameter α and within a typical range of bulk Lorentz factors, Γmin ≤ Γ ≤ Γmax, where Γmin = 1.034 and Γmax = 5.

Fig. 10. Fraction of the luminosity of protons transfered to neutrons for scenarios comprising different degrees of collimation of the jet.

In Fig. 11, we show the integrated spectrum of relativistic neutrons for models comprising different degrees of collimation of the jet. Again, all the spectra are broken power laws with similar spectral indices independently of the collimation regime. Moreover, the spectral indices agree well with those obtained in all previous scenarios. The energy at which the spectra break display smaller changes in this case, shifting to higher values as the jet is more collimated. This is consistent with our previous interpretation because more collimated jets show smaller adiabatic losses and, therefore, higher $E_{\text{p}}^{\text{max}}$ at the same z.

Fig. 11. Neutron power spectra in the jet frame, for different values of the collimation parameter α and Emin = 95.4mc2. The remaining parameters are the same as those in Table 1.

From the previous analysis, we can conclude that slow and collimated jets, with a compact region of acceleration close to the base and with negligible convection, may be very efficient sources of relativistic neutrons and, therefore, of cosmic rays.

3.3. Photohadronic interactions

The interaction of relativistic protons with radiation fields also provides a mechanism to produce relativistic neutrons via some of its branching channels (as discussed in Sect. 2). In this section, we assess the contribution of this interaction. We consider all the relevant radiation fields in a MQ, both internal and external to the jet. Given that in the jet frame the relativistic protons have an isotropic momentum distribution, we assume that the photomeson interactions are also isotropic, regardless of the direction of the target photons. In the first case, we adopt the synchrotron photon field radiated from both relativistic electrons and protons, in the local approximation of Ghisellini et al. (1985). In this case, the photon density can be computed as:

$$\begin{array}{c}\hfill n_{\mathrm{sy}}^{(\mathrm{i})}=\frac{q_{\mathrm{sy}}^{(\mathrm{i})}r(z)}{c{\epsilon }_{\mathrm{ph}}},\end{array}$$(14)

where i = e, p, for electrons and protons, respectively, ε_ph is the photon energy, and $q_{\mathrm{sy}}^{(\mathrm{i})}$ is the synchrotron emissivity ($[q_{\mathrm{sy}}^{(\mathrm{i})}]={\mathrm{s}}^{-1}{\mathrm{cm}}^{-3}$)¹. Regarding the external fields, we consider those of the companion star, the accretion disc, and the corona surrounding the black hole. For the stellar radiation, we adopt a black body with effective temperature T_⋆ ≈ 30 000 K in the case of the fiducial model (representing that of the O9.7 companion star of the Cygnus X-1 system). The photons emitted by the corona can be described with a power law of the form (e.g., Vieyro et al. 2012):

$$\begin{array}{c}\hfill n_{\mathrm{c}}({\epsilon }_{\mathrm{ph}})=Q_{\mathrm{c}}{\epsilon }_{\mathrm{ph}}^{-\beta }exp\left(\frac{{\epsilon }_{\mathrm{ph}}}{{\epsilon }_{\mathrm{ph}}^{\max }}\right),\end{array}$$(15)

where ${\epsilon }_{\mathrm{ph}}^{\max }\approx 150$ KeV. The field is normalised to the total power of the corona, assumed to be 1% of the Eddington luminosity of the black hole. For the black-hole mass, we adopt a value of M_BH ≈ 21 M_⊙, which corresponds to that of the Cygnus X-1 system (Miller-Jones et al. 2021). Finally, the disc emission is computed adopting a geometrically thin, optically thick accretion disc. The temperature profile is taken from (Romero & Vila 2014):

$$\begin{array}{c}\hfill T(R)={\left[\frac{3GM\dot{M}}{8\pi {\sigma }_{\mathrm{SB}}{R}^3}(1-\sqrt{\frac{R_{\mathrm{in}}}{R}})\right]}^{\frac{1}{4}},\end{array}$$(16)

where R is the radial coordinate of the disc, σ_SB is the Stephan-Boltzmann constant, and Ṁ is the accretion rate. The disc extends from R_in = 0.9R_c to R_out = 100R_in, where R_c = 35R_g is the radius of the corona and R_g the gravitational radius.

Figure 12 shows the loss rates for photohadronic interactions between protons and the relevant photon fields, at the base of acceleration region. In Fig. 13, we present the total loss rate for these interactions together with the rest of the losses for protons at the base, logarithmic midpoint, and top of the acceleration region. It is clear that the rate of photohadronic interactions falls well below that of proton-proton collisions, except near the maximum energy of protons at the base of the jet, where disc photons provide a rate similar to that of p–p. Near $E_{\mathrm{p}}^{\max }$, however, the amount of power released in neutrons is negligible because synchrotron and adiabatic losses are dominant and the proton population carries a low fraction of the total energy of the population (taking into account typical values of the injection spectral index, 1.5 ≲ p ≲ 2.4). As long as we consider regions further above the base of the jet, photohadronic losses become less relevant since the maximum energy of protons decreases (see also Fig. 5). Therefore, the photohadronic contribution to neutron production is negligible in every case.

Fig. 12. Loss rates for photohadronic interactions between protons and the photon fields of the companion star (blue), the accretion disc (green), the corona (orange), and local synchrotron emission in the jet (red). Dotted lines show the loss rate for the pair-creation channel and solid lines that of the photomeson channel. The total loss rate of the interaction (including pair-creation and photomeson channels) is plotted in solid purple line.

Fig. 13. Proton cooling rates including photohadronic interactions at the base (left panel), logarithmic midpoint (middle panel), and the end of the acceleration region (right panel). For this analysis, we use the same parameters as in Fig. 2.

3.4. Cosmic-ray spectrum

In this section, we compute the spectra of CRs produced as consequence of neutron decay far from the jet. Given the neutron energies in our model, the fraction of them that decay inside the jet is negligible. We focus on the secondary protons product of neutron beta decay. Electrons and neutrinos are also emitted in the process but they carry a negligible fraction of the progenitor neutron energy (∼0.1%). We compute the distribution of neutron decay distances using a probability exponential law in distance as in Escobar et al. (2021), with a mean equal to the lifetime of the neutron times its velocity, both in the ISM reference frame.

We compute the propagation of protons until they reach the ISM, assuming that they undergo elastic scattering with magnetic plasma fluctuations of the companion wind. We adopt as the boundary between the system and the ISM the sphere where the wind and external medium pressures become equal. We use the formulae given by Gleeson & Webb (1978) and Strauss et al. (2011) for the same process in the solar wind, which in our case yields:

$$\begin{array}{c}\hfill \dot{\gamma }=-\frac{2}{3}\gamma {\beta }^2\frac{v_{\mathrm{w}}}{r},\end{array}$$(17)

for the energy-loss rate of CR particles. In the previous expression, γ is the Lorentz factor of the particles, v_w is the stellar wind velocity, and r is the radial coordinate (centred in the binary system in this case). These energy losses are the result of particles propagating diffusively in the cavity through scattering off magnetic waves in the plasma. Adopting a diffusive radial motion with Bohm diffusion coefficient², the Lorentz factor of a proton that escapes to the ISM, γ_f, that had started to propagate with energy γm_pc² at a distance r₀ from the system is:

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{f}}-\gamma +\frac{1}{2}ln\left(\frac{{\gamma }_{\mathrm{f}}-1}{{\gamma }_{\mathrm{f}}+1}\right)=-\frac{2}{9}\frac{v_{\mathrm{w}}}{r_0^2}\frac{3e{B}_{\star }{R}_{\star }^3}{m_{\mathrm{p}}{c}^3},\end{array}$$(18)

where R_⋆ and B_⋆ are the star radius and magnetic field at the stellar surface, respectively. This equation allows us to compute the CR spectrum at the system boundary, assuming the number of particles is conserved throughout the propagation process.

In Fig. 14, we show the normalised CR spectra for all the models discussed in Sects. 3.1 and 3.2. In all cases, at energies of E ≳ 10¹¹ eV, the spectra are essentially the same as those of the neutrons. In addition, CR protons carry essentially all the neutron luminosity produced by the MQ (cf. Figs. 4, 7, and 10). This is because protons suffer almost negligible losses while propagating through the wind plasma. The only consequence of the interaction is the development of a low-energy tail in the spectra, arising from protons cooled below the minimum energy of the parent population. We discuss the implications of these findings in the next section.

Fig. 14. Cosmic-ray spectra for all scenarios discussed in this work, comprising different extensions of acceleration region (top), convection regimes (middle), and degrees of jet collimation (bottom). All spectra are normalised to the total CR luminosity.

4. Discussion and conclusions

In this work, we explore MQs as potential sources of CRs, whose production is mediated by relativistic neutrons created in the jet. We improve upon the work of Escobar et al. (2021) by taking into account the degree of collimation of the jet, the convective transport of particles, and the different sizes of the acceleration and emission regions. We also include photohadronic interactions with photon fields as another channel that could contribute to produce relativistic neutrons. In this way, we have explored a fairly general scenario for CR emission by a MQ jet, compatible with current observations.

Escobar et al. (2021) have shown that the jet luminosity and Lorentz factor are the main properties affecting the neutron luminosity and, therefore, the CR power injected into the ISM. They estimated the efficiency of the jets as CR sources in the range ∼10⁻⁶ − 10⁻⁴. In the present work, we show that, in addition, the collimation and compactness of the neutron production region, together with its proximity to the jet base, have a major influence on the efficiency of CR production. In the cases explored in this work, the CR power can rise close to four orders of magnitude for a highly collimated jet with a compact neutron emission region, in comparison with a conic jet. This renders MQ jets as serious candidates for an alternative CR acceleration engine. We note that the jets might even develop recollimation shocks that could yield a roughly cylindrical flow. In these cases, the efficiency of cosmic-ray production may increase even more, since the density of bulk protons stays nearly constant along the jet and so does the proton-proton collision rate. It is worth emphasising that the branching ratio of the neutron production channel in p–p collisions adopted in this work (0.16) is a low, conservative value and that the rate of neutron and CR production might increase up to a factor of ∼6 if this branching ratio is increased to the level adopted in previous works (Sikora et al. 1989; Atoyan 1992a,b; Vila et al. 2014; Romero & Gutiérrez 2020).

We recall that Pepe et al. (2015) modelled the Cygnus X-1 jet with a conic shape, whereas in a recent model by Kantzas et al. (2021), the authors assumed a well-collimated flow with no adiabatic losses. In both works, the authors manage to explain the available data, showing that current observations are still unable to break all model degeneracies. Therefore, for a single source, our model predicts a range of CR luminosities depending on the assumptions about the geometry of the acceleration region. For the case of Cygnus X-1, Escobar et al. (2021) estimate a CR luminosity of L_CR ≈ 10³² erg s⁻¹ (including the bulk and wind protons as targets for p–p collisions and the presence of a counterjet). In the present work, we show that for the same energetics of the relativistic populations, if the jet is assumed to be highly collimated, this amount would increase to L_CR ≈ 10^{35 − 36} erg s⁻¹. On the other hand, an ultraluminous X-ray source (ULX) would provide ∼10^{37 − 38} erg s⁻¹ if the jets are collimated. The total energy released in the MQ lifetime, adopting a duty cycle of ∼10⁶ yr for a high-mass MQ jet, such as that of Cygnus X-1, would be ∼10^{48 − 49} erg in the collimated case, whereas a ULX would directly rival a typical SNR, injecting 10^{50 − 51} erg. Clearly, an improvement in both multiwavelength observations and theoretical jet models is required to resolve these degeneracies and to estimate the contribution of MQ jets to the CR population with a higher precision.

Regarding the proton CR spectrum (see Fig. 14), our scenarios produce a broken power law with spectral indices of ∼2.3 (low energy) and ∼4.4 (high energy). The former is remarkably independent of any of the features of the jet (cf. Figs. 6, 8, and 11) and can be traced to the injection index, p, given that the loss mechanisms that dominate over the whole jet hardly modify the spectral index of the CR outcome. The break marks the maximum neutron energy at the top of the emission region. CRs with energies below the break are the result of the decay of neutrons coming from the full volume of the jet, whereas those above the break come progressively from smaller regions. This shrinking of the emission region produces the steepening of the CR spectrum at high energies. Both the break and the high-energy index show a deep connection with the basic physics of our models, mainly the distribution of magnetic field strength in the jet. Slowly decaying fields along the jet axis increase the energy range in which the whole region of emission produces neutrons, thus moving the break to higher energies and producing a steeper CR spectrum with less power above the break.

The high-energy part of the spectrum is softer than that predicted in the standard SNR scenarios. A harder injection spectral index would produce CR spectra with indices closer to that observed. Moreover, if we compute the propagation of CRs at energies in the range ∼10 GeV − 10 TeV (where the spectral index is ∼2.3) while adopting an ISM diffusion coefficient of D(E)∝E^0.3 (e.g., Blasi 2013), the spectra would evolve to a power law with a spectral index of ∼2.6, which is very close to that observed (∼2.7, e.g., Berezinskii et al. 1990); this is an aspect that the current SNR paradigm fails to explain. To fully settle this question, an estimate of the total spectrum of CRs produced by MQs, including the whole population of these sources and propagation effects, is required.

Even though photohadronic interactions are negligible in comparison with other proton losses in MQ jets, neutron production through these channel merits investigation in other scenarios that may provide suitable conditions for the process to take place significantly. This might be the case of MQ coronae, as explored in previous works (Vila et al. 2014). These scenarios would not only contribute to the total CR luminosity of MQs, but would also increase their duty cycle, as they would allow CRs to be produced even while a jet is not active. In fact, MQs may undergo phases in which a jet is present for a time, while over the rest of the cycle, the system may undergo different accretion episodes whereby the disc and coronae structures are formed.

Apart from quantitatively establishing MQ jets as serious candidates to compete with SNRs as sources of CRs, our work opens the possibility of taking a step forward in a different direction. Our results show that the CR emission depends only on jet properties. External photon fields and winds (either from the accretion disc, corona, or companion star) providing targets for photohadronic and p–p interactions, respectively, play a minor role in neutron production. They do not affect the propagation of neutron decay products until they reach the ISM either, except for very small changes at the low energy end of the spectrum. This indicates that only jet parameters are needed to predict the CR injection by a MQ, while the type of companion (i.e. high or low mass) and details of the disc or corona are irrelevant. By quantitatively relating jet parameters such as the Lorentz factor or the jet luminosity with the CR luminosity and spectrum, we pave the way to the construction of models describing the contribution of the whole population of MQs to the Galaxy CRs.

Finally, the mechanism explored in this work provides a tool for investigating the impact of MQs in the reionisation and heating rates in the early Universe. Different works have shown that X-ray binaries (and therefore MQs) would be more numerous and more luminous at Cosmic Dawn (Mirabel et al. 2011; Kaaret et al. 2011; Brorby et al. 2014; Douna et al. 2015; Ponnada et al. 2019). Our work provides a way of estimating the contribution of these systems to the CR population, whose role in the modification of the temperature and ionisation state of the intergalactic medium at early epochs is still a matter of debate (Tueros et al. 2014; Leite et al. 2017; Douna et al. 2018). In this context, our result on the existence of a low-energy tail of CRs produced in the propagation of neutron products in their way to the ISM is interesting given the larger ionisation power of low-energy particles (Douna et al. 2018).

Acknowledgments

G.J.E. and L.J.P. acknowledge the grant PIP 2014/0265 of Argentine National Research Council (CONICET). G.J.E. also acknowledges financial support from the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017. G.E.R. is supported by the National Agency for Scientific and Technological Promotion (PICT 2017-0898 and PICT 2017-2865) and the Spanish Ministerio de Ciencia e Innovación (MICINN) under grant PID2019-105510GBC31 and through the Center of Excellence Mara de Maeztu 2020-2023 award to the ICCUB (CEX2019-000918-M).

References

All Tables

All Figures

	Fig. 1. Schematic picture of the jet (not to scale), showing a jet launched at a distance z0 from the compact object. The acceleration region, shaded in red, is defined between za and zt, while the emission region spans over the range za ≤ z ≤ zmax.
In the text

	Fig. 2. Proton cooling rates for the case zt = zmax = 1015 cm, at the base (left panel), logarithmic midpoint (middle panel), and the end of the acceleration region (right panel).
In the text

	Fig. 3. Ratio of neutron emissivity to proton volumetric injection rate (top) and neutron luminosity distribution normalised with total luminosity Ln (bottom).
In the text

	Fig. 4. Fraction of the luminosity of relativistic protons transferred to neutrons, as a function of the size of the acceleration region, defined by the position of its top zt. Here, Lp stands for the total luminosity in the proton relativistic population.
In the text

	Fig. 5. Maximum energy of protons as a function of distance z along the jet for models comprising three different values of the acceleration region top zt. The top (bottom) panel corresponds to a magnetic index m = 1.9 (1.5).
In the text

Fig. 6. Integrated neutron power spectra in the jet frame, for different values of the top of the acceleration (and in this case also emission) region, zt = zmax. Blue stars indicate the value of Pn evaluated at the maximum neutron energy at zt = zmax, i.e. the minimum value of $E_{\mathrm{n}}^{\max }(=0.5E_{\mathrm{p}}^{\max })$ along the jet for each model (cf. Fig. 5). In all the cases Emin = 5mc2. The top (bottom) panel corresponds to a magnetic index m = 1.9 (1.5). The remaining parameters are those of Table 1.

In the text

	Fig. 7. Fraction of the luminosity of relativistic protons transferred to neutrons, as a function of the convective strength of the jet, defined by the parameter δ.
In the text

	Fig. 8. Integrated neutron power spectra in the jet frame, for scenarios with different convective strength, with Emin = 5mc2. The remaining parameters are the same as those in Table 1.
In the text

	Fig. 9. Ratio of proton–proton to adiabatic-plus-escape average loss rates, f, for models with different collimation parameter α and within a typical range of bulk Lorentz factors, Γmin ≤ Γ ≤ Γmax, where Γmin = 1.034 and Γmax = 5.
In the text

	Fig. 10. Fraction of the luminosity of protons transfered to neutrons for scenarios comprising different degrees of collimation of the jet.
In the text

	Fig. 11. Neutron power spectra in the jet frame, for different values of the collimation parameter α and Emin = 95.4mc2. The remaining parameters are the same as those in Table 1.
In the text

Fig. 12. Loss rates for photohadronic interactions between protons and the photon fields of the companion star (blue), the accretion disc (green), the corona (orange), and local synchrotron emission in the jet (red). Dotted lines show the loss rate for the pair-creation channel and solid lines that of the photomeson channel. The total loss rate of the interaction (including pair-creation and photomeson channels) is plotted in solid purple line.

In the text

	Fig. 13. Proton cooling rates including photohadronic interactions at the base (left panel), logarithmic midpoint (middle panel), and the end of the acceleration region (right panel). For this analysis, we use the same parameters as in Fig. 2.
In the text

	Fig. 14. Cosmic-ray spectra for all scenarios discussed in this work, comprising different extensions of acceleration region (top), convection regimes (middle), and degrees of jet collimation (bottom). All spectra are normalised to the total CR luminosity.
In the text