Introduction.— The presence of cosmic rays (CRs) indicates that the present universe has a high-energy aspect. Galactic CRs are thought to be accelerated by supernova remnants (SNRs) in our Galaxy. The gamma-ray telescopes AGILE and Fermi showed that Galactic CRs with energies below \(1~{\rm TeV}\) are actually accelerated by Galactic SNRs [1], [2]. In the standard picture of Galactic CRs, protons and irons are accelerated to about \(1~{\rm PeV}\) and \(100~{\rm PeV}\), respectively. However, it is not yet understood theoretically and observationally to what energies SNRs can accelerate CRs. The TALE experiment shows that a non-negligible amount of CR protons are present even at \(100~{\rm PeV}\) [3]. Theoretically, a strong magnetic field amplification or some special SNRs are required to accelerate protons to \(1~{\rm PeV}\) [4]–[7]. Moreover, the proton acceleration to \(100~{\rm PeV}\) requires both the strong magnetic field and a special SNR [8].

In addition to Galactic SNRs, black holes in our Galaxy have also been considered as CR sources [9]–[11]. It has been estimated that protons can be accelerated to \(100~{\rm PeV}\) in the Galactic black hole binary system (so called, microquasar) [11]–[14]. In fact, very high energy (VHE) gamma rays with energies exceeding \(100~{\rm TeV}\) from some microquasars was recently reported by LHAASO and HAWC [15]–[17]. The spatial extent of the VHE gamma rays is several tens of pc, which is much larger than the size of a black hole binary system, but comparable to a hot bubble produced by the wind or jet from the black hole binary [18]. Some microquasars have different gamma-ray spectra, some of which are steeper than predicted by the standard shock acceleration model [19], [20]. Although individual interpretations for the VHE gamma rays from the microquasars SS 433 and V4641 Sgr are given by [15]–[17], a unified explanation for the five microquasars that emit VHE gamma rays has not yet been attempted.

Gamma-ray spectra from middle-aged SNRs have similar characteristics: spatially extended gamma-ray distribution, different spectra, and steep spectra. These can be explained in a unified way by CRs escaping from SNRs [21]. Gamma-ray spectra from escaping CRs are steeper than predicted by the standard shock acceleration model due to the energy-dependent diffusion [22], [23]. In addition, the finite size of the escaping CR source has to be considered in order to understand the different gamma-ray spectra of some SNRs in a unified way [21], [24].

In this work, we try to understand the origin of the VHE gamma rays from microquasars in a unified way by considering CRs escaping from microquasars. For middle-aged SNRs, the CR acceleration to the TeV-PeV scale has already completed, so that the injection of escaping CRs can be treated as an impulsive injection. On the other hand, as long as the black hole continues to accrete mass, the microquasar would continue to accelerate CRs. Therefore, escaping CRs would be continuously injected into the surroundings. First, we derive the distribution of escaping CRs that are continuously injected in the finite size region. Then, we show that VHE gamma-ray spectra from five microquasars can be explained in a unified way by CRs escaping from the microquasars.

Energy spectrum of escaping CRs in a gamma-ray emitting region.— We consider a spherically symmetric system for simplicity. A microquasar launches a pair of jets or wind, creating a low-density hot bubble with the size of \(78~{\rm pc} (L/10^{39}~{\rm erg s^{-1}})^{1/5}(n/10~{\rm cm^{-3}})^{-1/5}(t/10^6~{\rm yr})^{3/5}\) [18]. CRs accelerated in the microquasar system are initially trapped in the bubble, but eventually escape from there. Therefore, gamma rays are expected to be produced by pp interaction outside the low-density hot bubble. We set \(R_1\) and \(R_2\) as the bubble radius and outer radius of the emission region, that is, gamma rays are produced in \(R_1\leq r\leq R_2\). Although the bubble expands with time, \(R_1\) and \(R_2\) are assumed to be constant in this work.

After escaping from the hot bubble region, the escaping CRs diffusively propagate to the surroundings. We derive the distribution of the escaping CRs at distance \(r\) from the microquasar center and at time \(t\), \(f(t,r,p)\), where \(p\) is the CR momentum. The Green function of the diffusion equation is given by \[G(t-t',|{\boldsymbol{r}}-{\boldsymbol{r}'}|,p)=\frac{e^{-\left(\frac{|{\boldsymbol{r}}-{\boldsymbol{r}'}|^2}{R_{\rm d}(t-t',p)}\right)^2}}{\pi^{3/2}R_{\rm d}(t-t',p)^3} ~~, \label{eq:fpoint}\tag{1}\] where \[R_{\rm d}(t-t',p) = \sqrt{4D_{\rm xx}(p)(t-t')} \label{eq:rd}\tag{2}\] is the diffusion length and \(D_{\rm xx}(p)\) is the diffusion coefficient, which is assumed to depend only on the CR momentum in this work. We assume that CRs are continuously escaping from the bubble radius, \(R_1\). Then, the source term in the diffusion equation is given by \[q_{\rm s}(t,r,p)=\frac{Q_{\rm esc}(p)}{4\pi r^2} \delta(r-R_1) ~~,\] where \(Q_{\rm esc}(p)\) is the escaping CR spectrum injected per time. Then, the solution to the diffusion equation is given by \[f(t,r,p)=\frac{Q_{\rm esc}(p)}{4\pi^{3/2}R_1r} \int_0^t {\rm d}t' \frac{e^{-\left(\frac{r-R_1}{R_{\rm d}(t-t',p)}\right)^2}-e^{-\left(\frac{r+R_1}{R_{\rm d}(t-t',p)}\right)^2}}{R_{\rm d}(t-t',p)}~~. \label{eq:f}\tag{3}\] Since the above integrand is approximately \(1/R_{\rm d}\) for \(r\sim R_1 > R_{\rm d}\), but \(4R_1r\exp\{-(r/R_{\rm d})^2\}/R_{\rm d}^3\) for the others [21], the approximate solution to the diffusion equation is \[f(t,r,p)=\frac{Q_{\rm esc}(p)t}{2\pi^{3/2}R_1R_{\rm d}(t,p)r} \label{eq:af1}\tag{4}\] for \(r\sim R_1 > R_{\rm d}\), and \[f(t,r,p)=\frac{Q_{\rm esc}(p)t}{\pi R_{\rm d}(t,p)^2r} {\rm erfc}\left(\frac{r}{R_{\rm d}(t,p)}\right) \label{eq:af2}\tag{5}\] for the others, where \({\rm erfc}(x) = (2/\sqrt{\pi})\int_x^{\infty}e^{-y^2}{\rm d}y\) is the complementary error function.

The spectrum of escaping CRs in the emission region \(R_1\leq r\leq R_2\) is given by \[\begin{align} F(t,p)&=& \int_{R_1}^{R_2} f(t,r,p)4\pi r^2{\rm d}r \nonumber \\ &=&Q_{\rm esc}(p)t M(t,p) \label{eq:fp} \end{align}\tag{6}\] where \(M(t,p)\) is given by

\({\rm erf}(x)=(2/\sqrt{\pi})\int_0^{x}e^{-y^2}{\rm d}y\) is the error function. If all escaping CRs are in the emission region, \(F(p,t) = Q_{\rm esc}(p) t\). Therefore, the function, \(M(p,t)\), describes the spectral modification due to the diffusive propagation of escaping CRs, which depends only on \(R_1/R_{\rm d}(t,p)\) and \(R_2/R_1\).

For \(R_2/R_1\sim 1\), since the emission region is close to the escape radius, the finite size of the escape region affects the spectrum of escaping CRs in the emission region, \(F(t,p)\). For \(R_{\rm d}(t,p) < R_2-R_1\lesssim R_1\), since about half of escaping CRs are in the emission region, \(M(t,p)\) is approximately given by \[M(t,p)\approx 0.5~. \label{eq:am1}\tag{8}\] For \(R_2-R_1<R_{\rm d}(t,p)<R_1\), CRs can diffuse beyond the emission region. Since propagation from \(R_1\) can be described by one-dimensional diffusion and Equation (4 ) can be applied, \(M(t,p)\) is approximately given by \[M(t,p)\approx \frac{R_2^2 -R_1^2}{2R_1\sqrt{\pi D_{\rm xx}(p)t} }~~. \label{eq:am2}\tag{9}\] For \(R_1<R_{\rm d}(t,p)\), since propagation from \(R_1\) can be described by three-dimensional diffusion and Equation (5 ) can be applied, \(M(t,p)\) is approximately given by \[M(t,p)\approx \frac{R_2^2 -R_1^2}{2D_{\rm xx}(p)t}~~. \label{eq:am3}\tag{10}\] Therefore, the spectrum of escaping CRs in the emission region, \(F(t,p)\), has two spectral breaks at \(p_{2-1}\) and \(p_1\) for \(R_1\sim R_2\), where \(p_{2-1}\) and \(p_1\) are given by the conditions \(R_{\rm d}(t,p) = R_2-R_1\) and \(R_{\rm d}(t,p)=R_1\), respectively. For \(p>p_1\), the steady state is realized, so that the CR spectrum in the emission region does not depend on the age of the system.

For \(R_2/R_1\gg 1\), the source region can be regarded as the point source because the escape region is much smaller than the emission region. In this case, from Equation (5 ), \(M(t,p)\) is approximately given by \[\begin{align} M(t,p) = \left\{ \begin{array}{ll} 1& ({\rm for}~R_{\rm d}(t,p)<R_2) \\ \frac{R_2^2 -R_1^2}{2D_{\rm xx}(p)t} & ({\rm for}~R_2<R_{\rm d}(t,p)) \end{array} \right. ~~. \label{rnxfjupg} \end{align}\tag{11}\] Thus, the spectrum of escaping CRs in the emission region, \(F(t,p)\), has one spectral break at \(p_2\), where \(p_2\) is given by the condition \(R_{\rm d}(t,p) = R_2\).

If the momentum dependences of the source spectrum and diffusion coefficient are given by \(Q_{\rm esc}(p)\propto p^{-s_{\rm esc}}\) and \(D_{\rm xx}(p)\propto p^{\delta}\), for \(R_2/R_1\sim 1\), the spectral index of escaping CRs in the emission region is approximately given by \[\begin{align} \frac{d\log F(t,p)}{d\log p} \propto \left\{ \begin{array}{lll} -s_{\rm esc}& (p<p_{2-1}) \\ -(s_{\rm esc}+\delta/2) & (p_{2-1}<p<p_1) \\ -(s_{\rm esc}+\delta) & (p_1<p) \\ \end{array} \right. , \label{eq:af3} \end{align}\tag{12}\] and for \(R_2/R_1\gg 1\), \[\begin{align} \frac{d\log F(t,p)}{d\log p} \propto \left\{ \begin{array}{ll} -s_{\rm esc}& (p<p_2) \\ -(s_{\rm esc}+\delta) & (p_2<p) \end{array} \right. ~~. \label{eq:af4} \end{align}\tag{13}\] The spectral index in the hardest momentum region represents the true spectral index of the total escaping cosmic rays, which is not always the same as the spectral index in the acceleration region [21], [25].

Figure 1 shows \(M(p,t)\) for some cases, where the diffusion coefficient, \(D_{\rm xx}\), is assumed to be the Bohm diffusion in \(10~ {\rm \mu G}\), \(t=10^5\) and \(10^7~{\rm yr}, R_1=30~{\rm pc}\), and \(R_2 = 40, 50\), and \(60 ~{\rm pc}\). For the Bohm diffusion coefficient, the characteristic break energy is estimated by \[cp_{\rm b} =0.2~{\rm PeV}\left(\frac{B}{10~{\rm \mu G}}\right)\left(\frac{R}{30~{\rm pc}}\right)^2\left(\frac{t}{10^5~{\rm yr}}\right)^{-1}~, \label{eq:ebreak}\tag{14}\] where \(R\) is \(R_2-R_1, R_1\), or \(R_2\). As can be seen in Figure 1, there are two characteristic breaks for \(R_2/R_1 \sim 1\) (\(R_2=40\) and \(50 {\rm pc}\)), and one characteristic break for \(R_2/R_1>1\) (\(R_2=60~{\rm pc}\)). The break energies and spectral indeces are consistent with Equations (12 )-(14 ).

Note that we have not taken account of the cutoff in the low energy side. If the inner radius of the emission region is larger than the bubble radius, low energy CRs with a small diffusion coefficient cannot propagate to the emission region for a finite time. As a result, the gamma-ray spectrum has a cutoff structure in the low energy side [21]–[23].

Gamma-ray spectrum from cosmic rays escaping from microquasars.—

Gamma-ray spectra from escaping CRs are calculated by using Equations (2 ), (6 ) and (7 ) and the Naima python package [26], in which gamma rays are calculated by the cross section of pp interactions given in [27]. By comparing our model with the observational data provided by LAHHSO and HAWC, we try to find a unified picture that the extended VHE gamma rays from five microquasars are produced by CRs escaping from the microquasars, rather than finding the best fitting parameters for each of the microquasars. The gamma-ray flux depends on the CR luminosity, distance, and total mass in the emission region. However, it is very difficult to estimate these values accurately. In this work, we adjust the normalization of the gamma-ray flux to fit the data.

Figures 2-6 show the VHE gamma-ray spectra for V4641 Sgr, SS 433, GRS 1915+105, MAXI J1820+070, and Cygnus X-1, respectively. Our model is in broad agreement with observational data. We find a unified picture that all of these objects have the same parameter values, except for their ages. In the unified model, the total energy spectrum of escaping CRs is \(E^{-2}\), the diffusion coefficient is the Bohm diffusion coefficient in \(10~{\rm \mu G}\), \(R_1=30~{\rm pc}\) and \(R_2=60~{\rm pc}\). The ages of these objects are assumed to be \(t=10^5~{\rm yr}\) for V4641 Sgr, \(t=10^6~{\rm yr}\) for SS433, and \(t=10^7~{\rm yr}\) for GRS 1915+105, MAXI J1820+070, Cygnus X-1. Therefore, the age of the objects is responsible for the spectral differences between them in the unified model. This interpretation of the VHE gamma-rays from microquasars is not unique because the spectral structure depends only on \(R_1/R_{\rm d}(t,p)\) and \(R_2/R_{\rm d}(t,p)\). Thus, the radius or the diffusion coefficient may also be responsible for the spectral difference.

Interestingly, our unified model suggests that for all microquasars, the total energy spectra of escaping CRs are \(dN/dE\propto E^{-2}\) with no cutoff at the high energy end. Since this spectrum is harder than that in SNRs [25], microquasars have the potential to supply Galactic cosmic rays above the PeV scale.

Discussion.— To explain the hard spectrum observed in V4641 Sgr and SS 433, \(R_1/R_{\rm d}=R_1(4D_{\rm xx}t)^{-1/2}\) should be about \(0.5\) and \(0.15\) at energy of \(cp=1~{\rm PeV}\). Thus, the diffusion coefficient at energy of \(cp=1~{\rm PeV}\) should be \[D_{\rm xx} =2.7\times 10^{27} {\rm cm^{2} s^{-1}} \left(\frac{R_1}{30~{\rm pc}}\right)^2 \left(\frac{t}{10^5~{\rm yr}}\right)^{-1}\left(\frac{R_1/R_{\rm d}}{0.5}\right)^{-2}. \label{eq:dxx}\tag{15}\] This value is the same as the Bohm diffusion coefficient in \(10~{\rm \mu G}\) and smaller than the Galactic mean value expected from the CR boron-to-carbon ratio [16], [17]. The smaller diffusion coefficient around the CR sources is also expected from gamma-ray observations of SNRs [28], [29]. Escaping CRs are expected to amplify the magnetic field turbulence through the CR streaming instability, so that the diffusion coefficient could be small around CR sources [30], [31]. Compared to escaping CRs around SNRs, escaping CRs can keep the diffusive flux sufficiently large around microquars for a long time because of the continuous injection. The self-confinement of CRs around microquasars should be studied in detail elsewhere.

In addition to CR protons, microquasars would accelerate CR electrons. Then, the CR electrons could escape from the microquasar system, producing extended gamma rays through inverse Compton scattering [32]. The propagation of CR electrons is different from that of protons because CR electrons suffer radiative losses. If we observe gamma rays from both CR electrons and protons and distinguish spatially, the diffusion coefficient around CR sources can be estimated from gamma-ray observations. We look forward to further observations by LHAASO and HAWC, and the future new gamma-ray and VHE gamma-ray observations, CTA, ALPACA, SWGO.

In this work, we have not investigated where and how cosmic rays are accelerated in the microquasar system. Many acceleration mechanisms are expected to work: shock acceleration in jet and wind [19], [20], [33], turbulent acceleration in the hot bobble [34], [35], shear acceleration at the jet edge [36]–[38], and so on. However, we do not understand which process mainly accelerates CRs in the microquasar system. The detailed physics of escape and acceleration of CRs should be investigated simultaneously. Since the microquasar system can be regarded as a smaller version of the active galactic nuclei system, we could understand the acceleration and escape processes of extragalactic CRs by observing the microquasar system. In terms of the steady-state injection of CRs, CRs accelerated in pulsar wind nebulae and stellar wind may also share common physics. Therefore, understanding CRs from microquasars is useful for understanding the high-energy aspect of the present universe.

Summary.— HAWC and LHASSO reported spatially extended VHE gamma rays from five Galactic black binaries (microquasars). The large emission region suggests that the VHE gamma rays originate from CRs escaping from the microquasars. To understand the distribution of escaping CRs around a microquasar, we have solved the diffusion equation taking into account the finite size of the CR source and the continuous CR injection. The energy spectrum in the emission region is described by a broken power law spectrum with one or two spectral breaks even though the total spectrum of escaping CRs is a single power law spectrum. We have found the unified picture that all five microquasars have the same energy spectrum of the escaping CRs, \(dN/dE\propto E^{-2}\), the same diffusion coefficient, and the same emission region. The hard energy spectrum of escaping CRs without the high energy cutoff supports the idea that the origin of Galactic CRs above the PeV energy is Galactic microquasars.

References↩︎