© ESO 2021

1. Introduction

Microquasars (MQs) are X-ray binaries (XRBs) with relativistic jets. These systems display a phenomenology that resembles that of active galactic nuclei (AGNs) but on smaller scales (Mirabel & Rodríguez 1994). As jets inject large amounts of energy into the interstellar region, they are expected to influence the surrounding medium. Tetarenko et al. (2018, 2020) present the detection of several transition spectral lines, which supports the presence of an interaction region between MQ jets and the interstellar medium (ISM). However, the authors find that the interaction region lies at much smaller scales than postulated in previous works. Microquasar jets may also heat the ISM by injecting kinetic energy. This issue has been investigated by Fender et al. (2005), who conclude that the input of kinetic energy into the ISM from MQ jets has a minor, but non-negligible contribution with respect to supernovae.

In addition, MQ jets could develop shocks that accelerate particles, injecting cosmic rays (CRs) into the Galaxy. Currently, the most plausible scenario is that supernova remnants (SNRs) accelerate CRs up to high energies. However, supernova shocks are in general non-relativistic, in contrast with the relativistic shocks in MQ jets. Thus, the spectrum of CRs produced in jets might be harder (i.e. with a larger fraction of high-energy CRs) than that of SNRs, and merits investigation.

Fender et al. (2005) estimate that 5 − 10 per cent of the CR power of the Galaxy might be produced in the MQ population. Likewise, Heinz & Sunyaev (2002) model a mechanism in which a narrow band of the CR spectrum is produced in the terminal shocks of MQ jets. This sharp spectrum has a characteristic proton energy that depends on the bulk Lorentz factor of the jet. Heinz & Sunyaev (2002) also argue that a broader CR spectrum may arise from the collective effect of MQs with different Lorentz factors. In a recent study, Cooper et al. (2020) reinforce the hypothesis of MQ jets as potential sources of CRs, suggesting that, because the maximum energy of jet-accelerated CRs is relatively high compared to other CR sources, the former may contribute significantly in the region between the Galactic and extragalactic components of the CR spectrum, that is, beyond the knee and below the ankle of this spectrum.

To inject CRs in the Galaxy, the scenarios devised in the aforementioned works (Heinz & Sunyaev 2002; Fender et al. 2005; Cooper et al. 2020) require the presence of hadronic matter in MQ jets. This is supported by some theoretical models (e.g., Blandford & Payne 1982), in which a large-scale magnetic field launches an outflow of material from the accretion disc by a magneto-hydrodynamic mechanism. This leads to jets composed of a hot plasma of thermal electrons and hadrons, with a relativistic component. On the other hand, some MQs were found to have hadronic content. Migliari et al. (2002) made observations of the jet of SS 433, revealing the presence of iron emission lines. Another case is the binary system 4U1630-47, for which Díaz Trigo et al. (2013) reported the detection of emission lines arising from baryonic matter travelling in the jet. Entrainment of matter from the stellar wind can also contribute to loading the jet with baryons (Romero et al. 2003).

In this paper, we explore a different mechanism by which MQ jets may inject kinetic energy and CRs into the ISM. Unlike previous models, in which CRs are produced at the terminal shock of the jet, we propose that CRs may be injected directly into the ISM by relativistic neutrons escaping from the jet. These neutrons are produced by the interaction of protons accelerated in internal shocks with thermal ones. Neutrons escape freely from the system because they do not interact with the magnetic field confining the plasma. They decay far away from the jet, injecting relativistic protons and electrons into the ISM.

The population of CRs produced in our scenario may differ from that of previous models for several reasons. Firstly, escaping neutrons carry almost all the energy of their relativistic proton progenitors, whilst particles in the jet undergo adiabatic and radiation losses in their way to the termination shock. For the same reason, the CR energy distribution of our model may be different from that of particles escaping through the termination shock. Secondly, the injection of CRs proceeds at distances much greater than jet scales. Thirdly, CR injection by neutrons is roughly isotropic, because MQ jets have small Lorentz factors. Finally, electrons arising in neutron decay may be more energetic than those accelerated at the termination shocks. This is because the energy distribution of the neutron population is related to that of the relativistic proton population in the jet, which is subject to smaller losses and therefore attains higher maximum energies.

Relativistic neutron production has already been explored in the jets of AGNs (Kirk & Mastichiadis 1989; Sikora et al. 1989; Giovanoni & Kazanas 1990; Atoyan 1992a,b; Atoyan & Dermer 2003). The latter authors compute the photo-meson production of relativistic neutrons and the γ-ray production from interactions of these neutrons with photon fields. They show that the production of narrow beams of ultra high-energy neutrons with E ≳ 10¹⁷ eV is possible. The existence of this neutron component has not yet been investigated in MQ jets. Such a population of relativistic neutrons may not only be a source of CRs, but may also have effects on the γ-ray emission of the jet.

This paper is organised as follows. In Sect. 2 we introduce the jet model, describing the relevant physical processes that drive neutron production. In Sect. 3 we use a fiducial set of parameters to apply the jet model. In addition, we compute the γ-ray spectrum of a typical MQ jet populated with protons and neutrons, and assess the observability of the neutron component. In Sect. 4 we follow neutron escape until decay, determining the power injected in CRs and their energy distribution. Finally, in Sect. 5 we discuss our results and present our conclusions.

2. Jet model

Our jet model is based on that of Romero & Vila (2008). Throughout this paper we use a fiducial set of parameters (summarised in Table 1) taken from Model A of Pepe et al. (2015). Pepe et al. (2015) model the electromagnetic output of the microquasar Cygnus X-1 accounting for the contribution of the companion star, jet, accretion disc, and a hot corona, and obtain two sets of parameters (models A and B) for their jet model that represent the best fit to the spectral energy distribution (SED) of Cygnus X-1. The uncertainties of the fitted parameters are not given, but we discuss the effects of parameter variations on our results in the following section.

Figure 1 presents a schematic picture of a MQ, showing the jets launched from the vicinity of the compact object. We adopt a lepto-hadronic, conic jet model with semi-aperture angle θ_jet. The flow propagates along the z axis with a macroscopic (bulk) Lorentz factor Γ, and carries a total power L_jet. The jet is pervaded by a magnetic field B(z) = B₀ (z₀/z)^α, where α is the magnetic index, and z₀ and B₀ are the position and magnetic field of the jet-launching region, respectively. The latter is derived assuming equipartition between magnetic and kinetic energy at the jet base,

$$\begin{array}{c}\hfill \frac{B_0^2}{8\pi }=\frac{L_{\mathrm{jet}}}{\pi r_0^2{v}_{\mathrm{jet}}},\end{array}$$(1)

Fig. 1. Top: schematic picture of a microquasar (not to scale). The main parameters used in the model are indicated. Bottom: discretisation mesh scheme defining a set of n logarithmically spaced points zi.

where r₀ is the jet radius at its base, and v_jet the jet bulk velocity. Hereafter we work in the jet frame, where the bulk is at rest.

We assume that internal shocks in the jet accelerate protons to relativistic energies in a region confined between z_min and z_max, and that cooling of these particles takes place locally. A fraction q_rel of the jet power is deposited in relativistic particles, with a proton-to-electron luminosity ratio a (i.e. a = L_p/L_e),

$$\begin{array}{c}\hfill q_{\mathrm{rel}}{L}_{\mathrm{jet}}=(1+a^{-1})L_{\mathrm{p}},\end{array}$$(2)

where L_p is the total power in the relativistic proton population. Injection of protons in the acceleration region is described by a source rate density

$$\begin{array}{c}\hfill Q(E_{\mathrm{p}},z)=Q_0{(E_{\mathrm{p}})}^{-p}exp[-E_{\mathrm{p}}/E_{\mathrm{p},\max }(z)],\end{array}$$(3)

where Q gives the number of protons injected per unit time and volume, and per proton energy interval (E_p, E_p + dE_p). The maximum attainable energy is E_p, max, p is the injection spectral index, and Q₀ is a normalisation constant obtained from

$$\begin{array}{c}\hfill L_{\mathrm{p}}={\displaystyle {\int }_V{\int }_{E_{\mathrm{p}}}Q(E_{\mathrm{p}},z)\mathrm{d}{E}_{\mathrm{p}}\mathrm{d}V,}\end{array}$$(4)

with V the volume of the acceleration region.

Neutrons are produced in inelastic collisions of relativistic protons with baryons in the bulk. We disregard the inverse process, that is, the production of protons by relativistic neutrons colliding with bulk protons. Below, we show that the rate of this process is negligible compared to the neutron escape rate. Photomeson production takes place beyond the threshold ϵ′≈145 MeV, where ϵ′ is the target photon energy in the frame of the relativistic proton. This corresponds to proton energies of ∼1 PeV for an X-ray photon field (e.g., that of the accretion disc or a corona). As we see below, contrary to AGN jets, these energy values are barely achieved by protons in a typical MQ jet. Other potential fields to be considered have lower-energy photons (e.g., synchrotron, stellar radiation) and therefore proton energies must rise above 1 PeV in order for the process to take place. Using a very similar model, Pepe et al. (2015, see their Fig. 5) found that the cooling rate for p–γ interaction against X-ray photons coming from a hot corona is several orders of magnitude lower than that for p–p collisions. Therefore, the production of neutrons through the former process is not taken into account.

In a stationary regime, the relativistic proton and neutron populations obey the following transport equations:

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

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

In these equations, b_p is the total energy loss rate of protons through all cooling processes (adiabatic plus radiative), and N_p(n) is the proton (neutron) spectral density. The term Q_n represents a source of neutrons due to proton–proton collisions, which is related to the proton sink term Λ through Q_n(E_n) = Λ(2E_n), as roughly half of the proton energy in a p–p collision is transferred to the neutron. The total neutron spectral power is then computed as

$$\begin{array}{c}\hfill P_{\mathrm{n}}(E_{\mathrm{n}})\equiv \frac{\mathrm{d}E}{\mathrm{d}{E}_{\mathrm{n}}\mathrm{d}t}={\displaystyle {\int }_V{E}_{\mathrm{n}}{Q}_{\mathrm{n}}(E_{\mathrm{n}})\mathrm{d}V,}\end{array}$$(7)

where [P_n] = s⁻¹. Neutrons experience neither adiabatic nor radiative cooling.

The last term of Eqs. (5) and (6) represents the escape of particles; following Romero & Vila (2008) we adopt $t_{\mathrm{esc}}^{-1}\sim {\mathit{v}}_{\mathrm{jet}}{(z_{\max }-z_{\min })}^{-1}$ as a characteristic escape rate for charged particles. Neutrons are not advected by the plasma, and therefore their escape rate is $t_{\mathrm{esc},\mathrm{n}}^{-1}\sim c/l$, with c the speed of light and l a characteristic size of the jet. A rather conservative lower limit is $t_{\mathrm{esc},\mathrm{n}}^{-1}>c{(z_{\max }-z_{\min })}^{-1}$, which arises from neutrons travelling through the whole acceleration region and escaping the jet through its head. Indeed, as neutron production is isotropic in the jet frame, most neutrons will escape by the side of the jet, increasing the rate by a factor of the order of ${\theta }_{\text{jet}}^{-1}$. As we show below, our key results do not depend on the exact value of $t_{\mathrm{esc},\mathrm{n}}^{-1}$.

Neutron decay is not considered at this stage because its rate is negligible with respect to escape. For our fiducial model z_max − z_min ≈ 2 × 10¹² cm, and we obtain $t_{\mathrm{esc},\mathrm{n}}^{-1}\gtrsim 1.5\times {10}^{-2}{\mathrm{s}}^{-1}$. The neutron decay rate is $t_{\mathrm{d}}^{-1}={\gamma }_{\mathrm{n}}^{-1}{\tau }_{\mathrm{d}}^{-1}$, where γ_n is the neutron Lorentz factor and ${\tau }_{\mathrm{d}}^{-1}=1.13\times {10}^{-3}\mathrm{s}$ is the inverse of the neutron lifetime. Therefore, neutron decay within the jet can be neglected in our model for any value of γ_n. This preserves the locality of the model, because neutron decay inside the jet would couple the populations of protons in different regions.

Radiative cooling processes considered for protons are synchrotron radiation due to the motion in the jet magnetic field and p–p interactions leading to γ-ray emission through pion decay. The synchrotron cooling rate $t_{\text{syn}}^{-1}$ is taken from Blumenthal & Gould (1970). For proton–proton inelastic scattering, we consider the two main branches

$$\begin{array}{cc}& p+p⟶p+n+{\pi }^{+}+l{\pi }^0+k({\pi }^{+}+{\pi }^{-}),\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill & p+p⟶p+p+l{\pi }^0+k({\pi }^{+}+{\pi }^{-}),\hfill \end{array}$$(9)

where l and k are the neutral and charged pion multiplicities, respectively. The total cooling rate of the process is given by

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

where K_pp and σ_pp are the p–p inelasticity and cross-section, respectively, v_p is the velocity of relativistic protons relative to that of bulk ones (we adopt v_p ≈ c), and $n_{\mathrm{b}}\approx L_{\mathrm{jet}}(1-q_{\mathrm{rel}}){[m_{\mathrm{p}}{c}^2{\mathit{v}}_{\mathrm{jet}}\pi {tan}^2{\theta }_{\mathrm{jet}}{z}^2]}^{-1}$ the bulk proton density in the reference frame of the jet (cf. Vila et al. 2012). We compute K_pp and σ_pp following Kafexhiu et al. (2014). A parameterisation of the inclusive cross section for the channel of Eq. (8) at low energies is also given by these authors. For consistency with that parameterisation and with measurements of the inclusive cross-section at higher energies (e.g., Engler et al. 1975; Flauger & Mönnig 1976; Adare et al. 2013; Adriani et al. 2018), we assume a probability of 0.16 as a low, conservative value for the production of a neutron in a proton–proton collision, that is, $\mathrm{\Lambda }(E_{\text{p}})=0.16N_{\text{p}}{t}_{\text{pp}}^{-1}$. The proton energy loss rate given by the interaction channel (9) is therefore

$$\begin{array}{c}\hfill b_{\mathrm{pp}}:=\frac{\mathrm{d}{E}_{\mathrm{p}}}{\mathrm{d}t}=-0.84{\sigma }_{\mathrm{pp}}{n}_{\mathrm{b}}{v}_{\mathrm{p}}{K}_{\mathrm{pp}}{E}_{\mathrm{p}}.\end{array}$$(11)

The neutron–proton collision rate is $t_{\mathrm{np}}^{-1}\sim t_{\mathrm{pp}}^{-1}$; therefore, neutrons will escape without interacting with bulk protons if $t_{\mathrm{pp}}^{-1}\ll t_{\mathrm{esc},\mathrm{n}}^{-1}$. As the opening angle is small, most neutrons will escape through the side of the jet, for which $t_{\mathrm{esc},\mathrm{n}}^{-1}\approx c/r(z)\approx 3\times {10}^3(z_{\min }/z){\mathrm{s}}^{-1}$, and the condition for escaping without interacting is fulfilled all along the region of interest.

Protons are also cooled through adiabatic losses, at a rate given by

$$\begin{array}{c}\hfill t_{\mathrm{ad}}^{-1}=\frac{2}{3}\frac{v_{\mathrm{jet}}}{z}.\end{array}$$(12)

Densities, cooling, source, and sink terms in the transport equations depend on particle energies and z. As there are no explicit spatial derivatives, Eqs. (5) and (6) become a set of coupled ordinary differential equations at each point z along the jet axis. To solve them, we discretise the functions in a mesh along the z axis in the region of interest. For each point in the mesh, we solve the system of equations numerically via the Picard method, using the one-zone approximation (i.e. assuming that what happens in one region has no effect on any of the others) and explicit differences to compute energy derivatives. Because the jet extends over several orders of magnitude in z, we adopt a logarithmic mesh (see Fig. 1). We use standard quadrature methods to compute the volume integrals required to obtain the properties of the whole acceleration region.

The boundary condition for solving the transport equations is N_p(E_p, max) = 0. The maximum energy of protons is obtained from the condition $t_{\mathrm{acc}}^{-1}(E_{\mathrm{p},\max })=t_{\mathrm{loss}}^{-1}(E_{\mathrm{p},\max })$, where the total loss rate is $t_{\mathrm{loss}}^{-1}=t_{\mathrm{ad}}^{-1}+t_{\mathrm{sync}}^{-1}+t_{\mathrm{pp}}^{-1}$, and $t_{\mathrm{acc}}^{-1}$ is the proton acceleration rate. We assume that a diffusive shock mechanism operates to accelerate charged protons at a rate of

$$\begin{array}{c}\hfill t_{\mathrm{acc}}^{-1}=\frac{\eta ecB}{E_i},\end{array}$$(13)

where e is the elementary charge, and 0 ≤ η ≤ 1 an efficiency parameter.

3. Neutron production and escape

3.1. The jet of Cygnus X-1

Figure 2 shows the proton acceleration and cooling rates for our fiducial model, at the base, middle, and top of the acceleration region, together with the neutron and proton densities at the same places. At the base, protons reach maximum energies of about 10¹⁵ eV, which is also an upper limit for the energies of neutrons, because the latter are about half of the former. At higher z, this value decreases because the acceleration rate, governed by the magnetic field, varies as z^−1.9, whereas the total loss rate changes as ∼z⁻¹. Adiabatic losses dominate along the whole jet. We therefore expect that the most energetic neutrons are mainly produced in regions near the base of the jet. On the other hand, we observe that the proton density is at least two orders of magnitude higher than the neutron density. The difference increases with the distance to the base of the jet. This is due to the fact that the density of target protons decreases with z. Thus, the proton spectrum is roughly the same as it would be without considering neutron production (we recall that the latter process is the unique sink for the proton population besides the escape), and the same applies to the gamma-ray SED produced by hadrons. This implies that, at least with the sensitivity of present instruments, there is no possibility of detecting the hadronic nature of jets by any signature produced by neutrons in their SEDs.

Fig. 2. Energetics of hadrons at the base (zmin, left panels), logarithmic midpoint (zmid, middle panels), and top (ztop, right panels) of the acceleration region for the Cygnus X-1 model. Top: cooling and acceleration rates for protons. The plots show the loss rates for proton–proton (green, solid line), synchrotron (yellow, dotted line), escape (magenta, +symbols), and adiabatic losses (light blue, dashed line), the total loss rate (black triangles), and the acceleration rate (red, dash-dotted line). Bottom: proton (black, solid line) and neutron (grey, dashed line) densities for the same regions as the top plots. In all cases, adiabatic losses are dominant and the proton population is many orders of magnitude denser than the neutron one.

The steady-state neutron density depends on the escape rate (see Eq. (6)), and represents the population of neutrons in transit before escaping. However, the injection rate of neutrons in the ISM is independent of the escape regime, because it is determined by the neutron production rate alone. Assuming that the neutron distribution is isotropic in the reference frame of the jet, the power injected in neutrons into the ISM, in the reference frame of the latter, will depend on the z-axis direction cosine μ in our model. To compute the spectral power of the neutron population injected into the ISM, we transform Eq. (7) according to

$$\begin{array}{c}\hfill P_{\mathrm{n}}^{\prime }(E_{\mathrm{n}}^{\prime })=\frac{1}{2}{\displaystyle {\int }_{-1}^1{P}_{\mathrm{n}}(E_{\mathrm{n}})\frac{\mathrm{d}\mathrm{\mu }}{\mathrm{d}{\mathrm{\mu }}^{\prime }}\frac{\mathrm{d}{E}_{\mathrm{n}}}{\mathrm{d}{E}_{\mathrm{n}}^{\prime }}\mathrm{d}{\mathrm{\mu }}^{\prime },}\end{array}$$(14)

where

$$\begin{array}{cc}& \frac{\mathrm{d}\mathrm{\mu }}{\mathrm{d}{\mathrm{\mu }}^{\prime }}=\frac{{(1+{\beta }_{\mathrm{jet}}{\mathrm{\mu }}^{\prime })}^2}{1-{\beta }_{\mathrm{jet}}^2},\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill & \frac{\mathrm{d}{E}_{\mathrm{n}}}{\mathrm{d}{E}_{\mathrm{n}}^{\prime }}=\mathrm{\Gamma }(1-{\beta }_{\mathrm{jet}}{\mathrm{\mu }}^{\prime }),\hfill \end{array}$$(16)

where β_jet = v_jet/c. In Eq. (14), the factor 1/2 comes from a previous integration in azimuth, and non-primed quantities refer to the rest frame of the jet. For simplicity, hereafter we use primed quantities for variables measured in the reference frame of the ISM.

Figure 3 shows the spectral power of the neutron population that escapes from the jet for the Cygnus X-1 case, which is computed in the reference frame of the ISM. The contributions of different regions of the jet are shown in the same way as in Fig. 2. The population produced at the base of the jet contributes to the high-energy spectral region, while that coming from higher zones dominates the low-energy region. The minimum energy is the same in all cases (E_n, min ≈ 0.5E_p, min) because that of protons is an input parameter, which for the fiducial model is the value determined by Pepe et al. (2015) via SED fitting. Regarding the total production, most of the power is injected in low-energy neutrons. The population presents a spectral index ($\stackrel{\sim }{p}$) of ≈3 at neutron energies of E_n ≈ 10^{11 − 12} eV, and steepens at higher energies, where the values of the spectral index shift up to $\stackrel{\sim }{p}\approx 3.5$ at E_n ≈ 10^{14 − 15} eV. The slope of the curve depends on the proton loss process that dominates at each energy. We note from Fig. 2 that, at each region of the jet, the relative contributions between adiabatic, proton–proton, and escape losses are modified. For the same reason, the maximum proton energies are also different in each region. Both effects contribute to a variation of the spectral index of the proton population, which is reflected in the spectral index of the neutron population, as shown in Fig. 3. The complete neutron population carries a total power of ≈3.3 × 10³¹ erg s⁻¹, which is of the order of 10⁻⁷L_jet.

Fig. 3. Spectral power of neutrons produced in collisions between relativistic protons with bulk ones computed in the ISM reference frame for the Cygnus X-1 model. We show the contribution of base, middle, and top regions in the jet (dashed, dash-dotted and dotted lines, respectively). The black solid line accounts for the total power. Violet segments represent power-law distributions with spectral index $\stackrel{\sim }{p}$.

The wind of the companion star may penetrate the jet and mix with its matter, thereby increasing its density and enhancing proton–proton interactions (Romero et al. 2003). The stellar-wind proton number density is given by

$$\begin{array}{c}\hfill n_{\mathrm{w}}=\frac{\dot{M}}{4\pi r^2{v}_{\mathrm{w}}{m}_{\mathrm{H}}},\end{array}$$(17)

where Ṁ is the mass-loss rate, v_w is the velocity of the wind, and m_H is the mass of the hydrogen atom. The standard velocity profile for a line-driven wind is given by (Lamers & Cassinelli 1999)

$$\begin{array}{c}\hfill v_{\mathrm{w}}=v_{\infty }{(1-\frac{R_{\star }}{r})}^{\delta },\end{array}$$(18)

where v_∞ is the terminal velocity, R_⋆ the radius of the star, and 0.5 ≤ δ ≤ 1.

For our fiducial model we adopt R_⋆ ≈ 20 R_⊙, v_∞ ≈ 2 100 km s⁻¹ and Ṁ ≈ 3 × 10⁻⁶ M_⊙ yr⁻¹ (e.g., Herrero et al. 1995; Yan et al. 2008), and the binary system separation a_* ≈ 3 × 10¹² cm (e.g., Iorio 2008). Assuming that all the material of the wind mixes with the jet, the bulk-to-wind density ratio is ≈8 × 10⁻²[1 + (a_⋆/z)²]. Thus, at the base of the acceleration region the bulk density overcomes that of the wind (n_p/n_w ≈ 10⁶), while near the top of the region the wind density becomes significant (n_p/n_w ≈ 10⁻¹). The contribution to the total neutron power is ≈6 × 10³¹ erg s⁻¹, and is roughly the same for any value of δ in the given range. This contribution is twice that of the neutron production with bulk protons as targets. We note that this wind scenario is extreme in the sense that the star is very close to the jet, and has a very high mass-loss rate. Varying the companion properties would then decrease the contribution of the wind material to neutron production. To achieve higher neutron luminosities, we therefore explore scenarios with different jet parameters.

3.2. Other jet scenarios

We have shown that, although a significative number of neutrons are indeed produced in a typical MQ jet, the energy carried by them to the ISM is small in the considered case. In this section we perform a variation of the main parameters of the jet model: the bulk Lorentz factor Γ, the efficiency of the acceleration η, the spectral index p, and the jet luminosity L_jet. The rest of the parameters remain those of the Cygnus X-1 model. In particular, we fix the value of the magnetic index (α = 1.9). The magnetic field is expected to have a dominant poloidal component near the jet base that becomes toroidal towards higher distances. This would result in a variation of the magnetic index from α ≈ 2 to α ≈ 1. However, as we see from Fig. 2, synchrotron losses dominate close to maximum energies and just near the jet base. In this region, the magnetic field is near the equipartition value, and therefore synchrotron may play a role in limiting the maximum neutron energy at the jet base only for jets with large kinetic luminosities. Beyond z ∼ z_mid, the synchrotron rate is negligible regardless of the value of the magnetic index. Figure 4 shows how the energetics of the hadron populations and the neutron production are modified in different scenarios. As we see, a more effective acceleration shifts the maximum energy of protons by ≳1 order of magnitude. On the other hand, lower bulk Lorentz factors lead to an increase in the proton–proton interaction rate, and hence the neutron production rate. For this parameter, we used alternative values, namely those of the MQ with the lowest jet velocity measured (Γ = 1.034, for SS 433; Chaty 2007) and the value adopted by Heinz & Sunyaev (2002, Γ = 5). Finally, we observe that increasing the jet luminosity increases both the relativistic and thermal proton densities, in turn increasing the neutron production rate.

Fig. 4. Left: total loss-rate (black solid line) and acceleration rates (red lines) for different values of the efficiency parameter η: the Cygnus X-1 model (solid line, η = 6 × 10−4), η = 10−2 (dashed line) and η = 0.1 (dotted line). Middle: acceleration rate (red solid line, η = 6 × 10−4), total loss rate (black lines), and proton–proton loss-rate (green lines) for different bulk Lorentz factors: that of Cygnus X-1 (solid line, Γ = 1.25), Γ = 1.034 (dashed line), and Γ = 5 (dotted line). Right: acceleration rate (red solid line) for the case of Cygnus X-1, and total loss rate (black lines) and proton–proton loss-rate (green lines) for different jet luminosities, i.e. Ljet: that of Cygnus X-1 (solid line, Ljet = 1038 erg s−1), Ljet = 1039 erg s−1 (dashed line), and Ljet = 1040 erg s−1 (dotted line). All panels refer to the jet base.

We show the spectral power of neutrons for eight models in Fig. 5. The variation of microscopic parameters such as η and p changes the hardness of the neutron population, increasing it as the proton injection becomes harder or the acceleration more efficient. These parameters produce minor variations in the total neutron power. The minimum energy of the relativistic protons is a parameter given in the rest frame of the jet; it changes the way in which the total energy input is distributed and is related to the value at which the neutron population peaks in the ISM frame. However, it does not have an impact on the spectral index or the shape of the distribution in general. On the other hand, macroscopic parameters do not significantly modify the spectral index, but do affect the general energetics of the population. A decrease in the bulk Lorentz factor increases the total power, while approximately preserving the shape of the spectrum. Higher Lorentz factors do not lead to significant changes in the spectrum because the lower densities of target protons limit the energy loss by the p–p channel. The neutron spectral power is also highly dependent on the jet luminosity, because an increase in the latter increases both the proton population energy and collision rate. For the eight models explored, we computed the total power in relativistic neutrons, L_n. The results are summarised in Table 2.

Fig. 5. Spectral power of relativistic neutron population for different model parameters, in comparison to that of Cygnus X-1. In all cases, the spectra are computed in the ISM reference frame. Left: variation of macroscopic parameters Γ and Ljet. Right: variation of microscopic parameters η and p.

Figure 6 shows the neutron-to-jet-power ratios for a wider range of parameter values. The variation of the spectral index and the efficiency parameter have a mild impact on the total neutron power. For reasonable values of these parameters, the power ratio varies by less than an order of magnitude. On the other hand, we observe a greater effect when varying the Lorentz factor and luminosity of the jet, which produces changes in the neutron power of several orders of magnitude. We note a rapid increase as Γ → 1. This is due to the increase in the bulk proton density, which results in a higher p–p rate. Therefore, as Γ increases, the p–p rate decreases, and so does the neutron power. At Γ ≈ 3, the effect of the Lorentz boost overcomes that of lower neutron production rates, producing a slight increase in the neutron power. Another important result is that the power ratio increases almost linearly with the jet luminosity. In other words, luminous jets are more efficient in transferring energy to the neutron component. This arises because the density of bulk and relativistic protons are both proportional to the jet luminosity, rendering the total neutron power L_n quadratic in L_jet. A general result of this section is that our model predicts that slow, high-luminosity jets are the astrophysical systems in which energy feedback into the ISM by neutron transport may play an important role.

Fig. 6. Neutron-to-jet-power ratios for different parameter values. Black dots correspond to the parameters of the Cygnus X-1 model and blue dots to the rest of the models. We perform variations of one parameter in each case, with the remaining parameters fixed to the values of the Cygnus X-1 model. From left to right, the variable parameter is: injection index p, acceleration efficiency η, the bulk Lorentz factor of the jet Γ, and its luminosity Ljet.

4. Cosmic-ray production

4.1. Neutron decay

We consider beta decay of free neutrons, $n⟶p+e+{\overline{\nu }}_e$ (Fermi 1934, et seqq.). The decay distance r follows an exponential probability density function f(r; γ_n) with mean cγ_nτ_n, where γ_n is the neutron Lorentz factor and τ_n ≈ 881.5 s is the neutron lifetime (Wietfeldt 2018). The power deposited within r in secondary particles is given by

$$\begin{array}{c}\hfill P_{\mathrm{d}}^{\prime }(r^{\prime })=\frac{1}{2}{\displaystyle {\int }_0^{r^{\prime }}{\int }_0^{\infty }{\int }_{-1}^{1}P(E_{\mathrm{n}})\frac{\mathrm{d}\mathrm{\mu }}{\mathrm{d}{\mathrm{\mu }}^{\prime }}\frac{\mathrm{d}{E}_{\mathrm{n}}}{\mathrm{d}{E}_{\mathrm{n}}^{\prime }}f(r^{\prime };{\gamma }_{\mathrm{n}}^{\prime })\mathrm{d}{\mathrm{\mu }}^{\prime }\mathrm{d}{E}_{\mathrm{n}}^{\prime }\mathrm{d}{r}^{\prime },}\end{array}$$(19)

In Fig. 7 we show P_d as a function of distance. We observe that, in our Cygnus X-1 model, most of the power is injected at distances ≳10¹⁵ cm and up to ∼10¹⁷ cm, values which are ∼10³ − 10⁵ times the binary system separation. This distance increases slightly, about half an order of magnitude, for models that produce harder neutron spectra.

Fig. 7. Power injected in secondary particles arising in neutron decay, Pd, normalised to Ln, up to a sphere of radius d′ centred at the MQ. Dashed horizontal lines indicate the first (blue), second (orange), and third quartiles of the total power. Four models presenting different behaviour are shown: Cygnus X-1 (black solid line), and models 4, 6, and 8 (red dashed, lilac dot-dashed, and brown dotted lines, respectively).

Secondary particles propagate through the matter, radiation, and magnetic fields surrounding the system. We consider the stellar wind of the companion as the main matter field in the decay region. This wind will expel the ISM matter and form a cavity of radius R_sys, which is given by the distance where the wind pressure is equal to that of the ISM (Fig. 8), beyond which we assume a typical ISM field. Thereby, we assume that particles that escape from the cavity become CRs. Given the high velocity and mass-loss rate of massive stars, it is expected that R_sys ≫ R_⋆, and therefore we can take v_w ≈ v_∞ for the velocity of the stellar wind in that region. Thus, the distance at which both pressures equilibrate is given by

$$\begin{array}{c}\hfill R_{\mathrm{sys}}={\left(\frac{\dot{M}{v}_{\infty }}{24\pi p_{\mathrm{ISM}}}\right)}^{1/2},\end{array}$$(20)

Fig. 8. Picture of the transport of particles within the stellar wind cavity (not to scale). The cavity is centred at the MQ. Neutrons propagate radially outwards until they decay into protons and electrons (neutrinos can be neglected for the purpose of this work). Charged particles follow a stochastic motion due to diffusion in the stellar wind plasma, losing energy until they reach the ISM and become CRs.

where p_ISM is the pressure of the ISM.

For the stellar wind of HDE 226 268 (the massive O9.7 star in Cygnus X-1 system) and a typical value of p_ISM ≈ 10⁻¹² dyn cm⁻², we obtain R_sys ≈ 2.3 × 10¹⁹ cm. Therefore the injected particles propagate inside the cavity formed by the stellar wind before escaping. The same applies for lower mass-loss rates, down to Ṁ ≈ 10⁻⁸ M_⊙ yr⁻¹, and for the whole range of wind velocities of massive stars (∼100 − 3000 km s⁻¹, e.g., Clark et al. 2012). Therefore, in high-mass MQs, neutron products would almost always have to travel some distance to reach the ISM, losing part of their energy.

The energy deposited in neutron-decay products may be carried away from the system by them, radiated through their interactions with magnetic, photon, and matter fields, or transferred to the medium by elastic interactions. Adopting typical values for stellar magnetic fields (B ≈ 100 G for the surface of a high-mass star), synchrotron losses are negligible in comparison to the adiabatic losses that particles suffer when propagating through the wind plasma. Relevant fields for proton–photon or electron inverse Compton interactions come from the binary system (companion star, accretion disc, jet, etc.). However, the collision rate is negligible in both cases because the encounter is produced at very small angles (≲0.001), as both colliding particles propagate outwards away from the system. Regarding proton–proton inelastic collisions, the mean free path is ≳10 pc for the wind-matter field. On the other hand, for typical values of magnetic field, the electron synchrotron cooling rate implies that their energy is radiated within typical distances ≳1 kpc, depending on the neutron decay distance. Thereby, the emission inside the cavity would be negligible. Thus, radiative losses of these particles are negligible while they propagate towards the ISM. Instead, before emerging as cosmic rays, part of their energy is lost while diffusing through the plasma.

4.2. Cosmic-ray spectra

To compute the losses of secondary particles in the stellar wind until they reach the edge of the cavity, we take the formulae for the same process in the solar wind, given by Gleeson & Webb (1978) and Strauss et al. (2011). In our case, the energy-loss rate can be written as

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

where γ is the Lorentz factor of the particle and v_w the stellar wind velocity at a given distance r from the binary system. These energy losses are the result of particles propagating diffusively in the cavity through scattering off magnetic waves in the plasma.

The radial motion equation of relativistic particles is given by $r=\sqrt{Dt}+r_0$, where r₀ is the injection –neutron decay– distance and D is the diffusion coefficient, for which we adopt the Bohm approximation. In terms of the model parameters, $D\approx Ec{r}^3{r}_{\star }^{-3}/(3e{B}_{\star })$, where B_⋆ and R_⋆ are the surface magnetic field and radius of the companion star, respectively. We use this relation to integrate Eq. (21) from the injection position r₀ to R_sys, yielding

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{F}}-\gamma (r_0)+\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}{\mathrm{m}{\mathrm{c}}^3},\end{array}$$(22)

where γ_F = γ(R_sys). The value of γ(r₀) depends on the decaying neutron energy, which is distributed among the created proton (99.9%) and electron (0.1%). Equation (22) gives the energy at which particles escape the system and emerge as cosmic rays in the ISM. Using this equation, we compute the cosmic-ray spectra assuming the number of particles is conserved for each population.

Figure 9 shows the final Lorentz factor of secondary particles injected at different distances. Particles coming from neutrons that decay near the source suffer more losses as they diffuse through a longer path. The proton population does not suffer significant losses for decay distances greater than ∼10¹⁴ cm. The effect of diffusion is greater for electrons, even at large decay distances. Figure 10 shows the spectra of CRs (protons and electrons) at r = R_sys for our fiducial model. The cosmic-ray proton spectrum is almost the same as the injected one, with a spectral index of p ∼ 3, but presents a tail at low energies due to diffusion. The electron spectrum flattens as a consequence of electrons suffering higher diffusion losses, and accumulating at lower energies in a tail similar to that of the proton spectrum. Both spectra show a maximum value around, and related to, the minimum Lorentz factor of the proton population. We recall that for our fiducial model, E_min = 95.4 mc². For other jet models, this value may vary down to ∼2 mc², and consequently, the maximum of the CR spectra will be located at lower energies.

Fig. 9. Lorentz factor at escape, γF, vs. Lorentz factor at injection, γ, for protons (top) and electrons (bottom) injected at different distances: d = 1014 cm (blue solid line), d = 1015 cm (dashed green line), d = 1016 cm (dashdotted brown line), and d = 1017 cm (dotted red line). For these results we adopted a typical surface value of B⋆ ≈ 100 G for a high-mass star.

Fig. 10. Cosmic-ray spectra produced by neutrons in MQs. Blue and red dashed lines represent injected electron and proton spectra, respectively. Black and brown solid lines represent the cosmic-ray spectra for the case of a companion star with B ≈ 100 G and B ≈ 3000 G, respectively. In the case of protons, all spectra are almost identical for energies E ≳ Emin.

The total CR spectrum carries essentially all the power deposited in the neutron population. Therefore, for any model explored in previous sections, the total cosmic-ray power is given by L_CR  ≈  L_n (see Table 2). Losses suffered by particles before escaping from the system depend on the wind velocity and injection distance. Therefore, other parameters will not produce significant variations for these energy losses, and the cosmic-ray spectra will modify analogously to the neutron spectra.

5. Discussion and conclusions

We introduced the relativistic neutron component in hadronic jet models through inelastic proton-proton collisions that produce these particles in situ. The density of the proton population overcomes that of the neutrons by a factor of ≳10² and neutron decay is negligible inside the jet. Therefore, the steady-state proton population is almost the same as that obtained without considering neutron production. The same is true for the jet SED, making the identification of the neutron component in MQs by its emission unattainable with present instruments.

Neutrons escape and decay far away from the jet, but remain inside the cavity carved out in the ISM by the MQ companion, as far as high-mass MQs are concerned. These neutrons inject secondary relativistic protons and electrons that propagate diffusively and finally escape from the cavity into the ISM, becoming CRs. These particles do not radiate significantly within the cavity, being undetectable through electromagnetic emission. They may effectively carry a small fraction of up to 10⁻⁴ of the jet power out into the ISM, depending mainly on the jet luminosity and bulk velocity. The distribution of this power among the proton and electron components is governed mainly by neutron decay physics, except by the small amount of power lost by electrons in diffusing out of the system.

Microquasars have been considered as potential CR sources in previous works (e.g., Heinz & Sunyaev 2002). These authors computed the CR output of cold protons and heavy ions that accelerate and escape through the terminal shock in the jet. The main feature of their CR spectra is a narrow shape around a typical Lorentz factor of at most a few times that of the jet bulk, which is usually Γ ∼ 3 − 10. Alternatively, our mechanism produces broad spectra peaked at half the minimum energy of relativistic protons in the jet, which may be more than one order of magnitude higher. The shape of the CR spectrum is similar for protons and electrons in our case; it presents an almost flat tail below the peak energy, and a steep decay with an index around −3 above. Briefly, our neutron-escape-based scenario can provide more energetic CRs because it drains energy directly from the relativistic proton population at the base of the jet, instead of taking that of cool particles emerging at its end.

It is important to recall that Heinz & Sunyaev (2002) assume that CRs emerging from the terminal shock of the jet are injected directly into the ISM. However, typical jets have lengths of up to ∼10¹⁵ cm, much smaller than the sizes of the cavities carved out by MQ companions in the ISM. Therefore, it is expected that the spectra of CRs exiting the jet through the terminal shock are modified by the wind. For particles with small Lorentz factors, such as those obtained by Heinz & Sunyaev (2002), and stars like the companion of Cygnus X-1, our results suggest that CRs would thermalise within the stellar wind if the terminal shock lies at ≲10¹⁴ cm. Therefore, the CR spectrum may be highly modulated by the wind, depending on the specific properties of the system. This issue merits a more thorough investigation.

Supernova remnants are at present considered the main sources of Galactic CRs. They inject ∼10⁵¹ erg into the ISM. However, the efficiency of the acceleration of CRs in SNR shocks remains under discussion. It is accepted that if ∼10% of the SNR energy is used in CR acceleration, the supernova population might explain the observed CR power in the Galaxy up to energies lower than that of the knee of the CR spectrum. Hovey et al. (2018) measured an upper limit of ∼7% for this efficiency, depending on the ionisation factor of the pre-shock gas, whereas Shimoda et al. (2015) argue that the CR acceleration efficiency may have been overestimated by 10 − 40%. Regarding our results, and assuming that a MQ jet like Cygnus X-1 might be active over a time of ∼10⁶ yr, it would inject only ∼10⁴⁵ erg into the ISM, far below the contribution of an individual SNR. Even slow jets, such as that of SS 433, would inject only 3 × 10⁴⁶ erg, still a low CR power. Only ultraluminous X-ray sources with the most powerful stellar jets can compete with SNRs by producing up to ∼10⁴⁹ erg, about 10% of the CR power of a SNR. More optimistic scenarios than those mentioned could be obtained combining parameters that favour neutron production independently, but would hardly represent the typical MQ population. It is important to recall that our estimates rely on a conservative value for the branching ratio of the neutron production channel. Other authors adopt values that are higher by a factor of up to six (cf. Sikora et al. 1989; Atoyan 1992a,b; Vila et al. 2014; Romero & Gutiérrez 2020), which would increase the neutron power by a similar amount.

On the other hand, measurements of the CR spectrum observed at Earth suggest that it is steeper than that predicted by actual models of diffusive acceleration in SNR shocks (e.g., Blasi 2013). In particular, spectral indices softer than ∼2 (the canonical value in the standard theory of diffusive shock acceleration) are required to describe observations. Our results suggest that these values may be easily achieved by a neutron-escape-based mechanism. This is due to a steepening of the neutron spectrum, which is the result of a contribution of several spectra with different values of the maximum energy achieved at each region in the jet (compare with Fig. 3).

A recent lepto-hadronic model for Cygnus X-1 was proposed by Kantzas et al. (2021). The main differences with the model on which ours is based (Pepe et al. 2015) are that the acceleration region extends to higher distances along the jet axis, and that synchrotron photons are adopted as targets for proton–photon interactions. This latter interaction is another source of relativistic neutrons that could dominate over the proton–proton channel in some cases. According to our model, the neutrons produced at farther regions would contribute to the lower energies of the population and modify the spectrum accordingly (see Fig. 3).

A more precise estimate of the contribution of MQs to Galactic CRs should take into account population considerations, because of the different production rates, lifetimes, and duty cycles of both classes of objects, which depend on the properties of the parent stellar populations. An interesting by-product is that old stellar populations might contribute to Galactic CRs through low-mass MQs. In the case of a low-mass companion star, harder neutron spectra could be expected, because their slow winds produce small cavities, allowing neutrons to inject a significant amount of power directly outside. This may also be the case for low-metallicity stars, which produce weaker winds. In these cases, the products of neutron decay emerge directly as cosmic rays. An exploration of these issues requires the development of stellar evolution models that include the MQ phase, and a more complete census of Galactic MQs to determine their actual population. A rough computation has already been presented by Fender et al. (2005), who estimate the MQ contribution to the Galactic cosmic-ray population to be in the range of 5 − 10%. Similar studies may also shed light onto the population of CRs of star-forming galaxies, and the origin of their γ-ray emission (see e.g., Romero et al. 2018; Kornecki et al. 2020, and references therein).

Based on the increase in the XRB production rate and luminosity at low metallicities, Mirabel et al. (2011) proposed that these sources may have contributed to the ionisation and heating of the IGM in the early Universe through their X-ray emission. This argument has been extended to include the contribution of CRs from MQs (Tueros et al. 2014; Douna et al. 2018). The latter authors emphasise that the electron spectrum is a key ingredient, finding that MQs with soft electron spectra provide the largest ionisation powers. Our work therefore provides a mechanism to support the claims of Douna et al. (2018). Moreover, Sotomayor Checa & Romero (2019) present theoretical models in which Population III MQs produce CRs in the terminal regions of extremely powerful hadronic jets (L_jet ∼ 10⁴¹ erg s⁻¹). Our neutron-escape-based mechanism predicts that those systems would have a very strong CR emission of ∼10³⁸ − 10³⁹ erg s⁻¹, and therefore a large ionising and heating power. Taking into account the fact that theoretical models predict that Population III stars have weak winds as a consequence of their low metallicities, relativistic neutrons would decay directly in the ISM, without suffering diffusion losses. Therefore, the contribution of our mechanism to ionising CRs in the early Universe merits exploration. We will present our results on this issue in a forthcoming paper.

Acknowledgments

The authors acknowledge the anonymous referee for valuable comments that greatly improved the manuscript. G. J. E. and L. J. P. acknowledge support from project PIP 2014/0265 from Argentine CONICET.

References

All Tables

All Figures

	Fig. 1. Top: schematic picture of a microquasar (not to scale). The main parameters used in the model are indicated. Bottom: discretisation mesh scheme defining a set of n logarithmically spaced points zi.
In the text

Fig. 2. Energetics of hadrons at the base (zmin, left panels), logarithmic midpoint (zmid, middle panels), and top (ztop, right panels) of the acceleration region for the Cygnus X-1 model. Top: cooling and acceleration rates for protons. The plots show the loss rates for proton–proton (green, solid line), synchrotron (yellow, dotted line), escape (magenta, +symbols), and adiabatic losses (light blue, dashed line), the total loss rate (black triangles), and the acceleration rate (red, dash-dotted line). Bottom: proton (black, solid line) and neutron (grey, dashed line) densities for the same regions as the top plots. In all cases, adiabatic losses are dominant and the proton population is many orders of magnitude denser than the neutron one.

In the text

Fig. 3. Spectral power of neutrons produced in collisions between relativistic protons with bulk ones computed in the ISM reference frame for the Cygnus X-1 model. We show the contribution of base, middle, and top regions in the jet (dashed, dash-dotted and dotted lines, respectively). The black solid line accounts for the total power. Violet segments represent power-law distributions with spectral index $\stackrel{\sim }{p}$.

In the text

Fig. 4. Left: total loss-rate (black solid line) and acceleration rates (red lines) for different values of the efficiency parameter η: the Cygnus X-1 model (solid line, η = 6 × 10−4), η = 10−2 (dashed line) and η = 0.1 (dotted line). Middle: acceleration rate (red solid line, η = 6 × 10−4), total loss rate (black lines), and proton–proton loss-rate (green lines) for different bulk Lorentz factors: that of Cygnus X-1 (solid line, Γ = 1.25), Γ = 1.034 (dashed line), and Γ = 5 (dotted line). Right: acceleration rate (red solid line) for the case of Cygnus X-1, and total loss rate (black lines) and proton–proton loss-rate (green lines) for different jet luminosities, i.e. Ljet: that of Cygnus X-1 (solid line, Ljet = 1038 erg s−1), Ljet = 1039 erg s−1 (dashed line), and Ljet = 1040 erg s−1 (dotted line). All panels refer to the jet base.

In the text

	Fig. 5. Spectral power of relativistic neutron population for different model parameters, in comparison to that of Cygnus X-1. In all cases, the spectra are computed in the ISM reference frame. Left: variation of macroscopic parameters Γ and Ljet. Right: variation of microscopic parameters η and p.
In the text

Fig. 6. Neutron-to-jet-power ratios for different parameter values. Black dots correspond to the parameters of the Cygnus X-1 model and blue dots to the rest of the models. We perform variations of one parameter in each case, with the remaining parameters fixed to the values of the Cygnus X-1 model. From left to right, the variable parameter is: injection index p, acceleration efficiency η, the bulk Lorentz factor of the jet Γ, and its luminosity Ljet.

In the text

Fig. 7. Power injected in secondary particles arising in neutron decay, Pd, normalised to Ln, up to a sphere of radius d′ centred at the MQ. Dashed horizontal lines indicate the first (blue), second (orange), and third quartiles of the total power. Four models presenting different behaviour are shown: Cygnus X-1 (black solid line), and models 4, 6, and 8 (red dashed, lilac dot-dashed, and brown dotted lines, respectively).

In the text

Fig. 8. Picture of the transport of particles within the stellar wind cavity (not to scale). The cavity is centred at the MQ. Neutrons propagate radially outwards until they decay into protons and electrons (neutrinos can be neglected for the purpose of this work). Charged particles follow a stochastic motion due to diffusion in the stellar wind plasma, losing energy until they reach the ISM and become CRs.

In the text

	Fig. 9. Lorentz factor at escape, γF, vs. Lorentz factor at injection, γ, for protons (top) and electrons (bottom) injected at different distances: d = 1014 cm (blue solid line), d = 1015 cm (dashed green line), d = 1016 cm (dashdotted brown line), and d = 1017 cm (dotted red line). For these results we adopted a typical surface value of B⋆ ≈ 100 G for a high-mass star.
In the text

	Fig. 10. Cosmic-ray spectra produced by neutrons in MQs. Blue and red dashed lines represent injected electron and proton spectra, respectively. Black and brown solid lines represent the cosmic-ray spectra for the case of a companion star with B ≈ 100 G and B ≈ 3000 G, respectively. In the case of protons, all spectra are almost identical for energies E ≳ Emin.
In the text