Abstract
Microquasars are compact binary systems hosting collimated relativistic jets. They have long been proposed as cosmic-ray accelerators, probed via the gamma-ray emission produced by relativistic particles. However, the observational evidence is steadily increasing but limited: there are around 20 microquasars known to date, of which only 3 have so far been firmly detected in the GeV gamma-ray range, always in a flaring or special spectral state. Here we present Fermi-LAT observations of the region around the microquasar GRS 1915+105, which reveal the presence of previously unknown multi-GeV emission consistent with the position of the microquasar. No periodicity or variability is found, indicating a persistent source of gamma rays. The properties of the emission are consistent with a scenario in which protons accelerated in the jets interact with nearby gas and produce gamma rays. We find that if the jet has been operating at an average of 1% of the Eddington limit for 10% of the time that GRS 1915+105 spent in its current mass transfer state, the transfer of 10% of the available power to protons is enough to reach the ∼3 × 1049 erg required to explain the GeV signal. Therefore, our results support a scenario in which microquasars with low-mass stellar companions act as hadronic accelerators, strengthening the idea that microquasars as a class contribute to at least some fraction of the observed cosmic-ray flux.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
1. Introduction
Stellar binary systems hosting a compact object reveal themselves through the X-ray emission associated with accretion onto the compact object, earning them the name of X-ray binaries. Microquasars are a subclass of X-ray binaries in which a collimated relativistic outflow of plasma (jet) is formed as a consequence of the aforementioned accretion (I. F. Mirabel & L. F. Rodrìguez 1999; R. P. Fender et al. 2004).
The microquasar GRS 1915+105 was first detected as an X-ray source by the WATCH instrument on board the GRANAT observatory in 1992 (A. J. Castro-Tirado et al. 1992). Follow-up observations with the Very Large Array (I. F. Mirabel & L. F. Rodrìguez 1994; L. F. Rodrìguez & I. F. Mirabel 1999) and MERLIN (R. P. Fender et al. 1999) in the radio band revealed a highly variable counterpart with apparently superluminal two-sided ejections. This was the first observation of relativistic motions in an object inside our Galaxy and implied intrinsic velocities for the ejecta near the speed of light. These findings established the presence of a jet with velocity v ∼ 0.8c and angle to the line of sight θ ∼ 63° (R. P. Fender et al. 1999; L. F. Rodrìguez & I. F. Mirabel 1999), making the system a microquasar.
After some initial constraints on the distance by R. P. Fender et al. (1999) and C. Chapuis & S. Corbel (2004), recent parallax measurements result in a distance estimate of d = 9.4 ± 0.6 ± 0.8 kpc (M. J. Reid & J. C. A. Miller-Jones 2023), which we will adopt through this work.
The mass of the black hole has been subject to debate, with an initial claim of MBH = 14 M⊙ (J. Greiner et al. 2001a), which was later revised to lower values ranging between 10 M⊙ (D. Steeghs et al. 2013) and 12 M⊙ (M. J. Reid et al. 2014; M. J. Reid & J. C. A. Miller-Jones 2023). The corresponding Eddington luminosity is on the order of LEdd = 1.26 × 1038MBH erg s erg s−1. The spectral type of the donor star was identified to be a K-type star (J. Greiner et al. 2001b) with a mass of 0.5 M⊙ (D. Steeghs et al. 2013), making GRS 1915+105 a low-mass X-ray binary (LMXB). Its orbital period is 33.85 ± 0.16 days (D. Steeghs et al. 2013).
Since its detection more than 30 yr ago, GRS 1915+105 has been extensively monitored from radio to X-ray energies, revealing a highly variable source with distinct behaviors across wavelengths. The emission in the radio band is characterized by quasi-periodic oscillations (G. G. Pooley & R. P. Fender 1997; R. P. Fender & G. G. Pooley 2000), seen also in the infrared band (R. P. Fender & G. G. Pooley 1998). The source flares episodically in both radio and X-rays, with the same anticorrelation between these two bands that is seen for other microquasars. Most notably, GRS 1915+105 has been a consistently bright X-ray source since its discovery, implying intrinsic X-ray luminosities in the range of 0.1LEdd–1LEdd, until a sudden decrease in 2018 June. An initial decrease brought the flux down by an order of magnitude, followed by a continuous decline of the X-ray emission while flaring activity increased at other wavelengths. This behavior has been attributed to either an obscuration (M. Balakrishnan et al. 2021; S. E. Motta et al. 2021) or an intrinsic change into a hard state (R. P. Fender et al. 2004; K. I. I. Koljonen & T. Hovatta 2021).
Gamma-ray emission from microquasars has long been predicted (e.g., A. M. Atoyan & F. A. Aharonian 1999; F. A. Aharonian 2004; V. Bosch-Ramon & D. Khangulyan 2009), in both the high-energy (HE; GeV-scale) and very high energy (VHE; TeV-scale) bands of the spectrum, usually invoking particle acceleration inside the jets. Only two objects, both high-mass X-ray binaries, have so far been firmly detected as GeV sources: Cyg X-3 (Fermi LAT Collaboration et al. 2009; M. Tavani et al. 2009; S. Corbel et al. 2012), where orbital gamma-ray modulation is seen, and Cyg X-1 (R. Zanin et al. 2016; A. A. Zdziarski et al. 2017), where gamma-ray emission is detected during its hard X-ray state. Additionally, a source has been detected near SS 433 (J. Li et al. 2020), showing flux variability correlated with the jet precession period. At higher energies only two microquasars have been firmly detected in the VHE band, both as persistent sources: SS 433 (A. U. Abeysekara et al. 2018; Z. Cao et al. 2024; H. E. S. S. Collaboration et al. 2024) and the intermediate-mass X-ray binary V4641 Sgr (R. Alfaro et al. 2024). No significant VHE emission has been reported from any of the other HE sources (M. L. Ahnen et al. 2017; A. Albert et al. 2021).
In comparison, the gamma-ray detection of LMXBs remains elusive despite intensive searches (particularly during the X-ray flaring states of V404 Cygni; A. Loh et al. 2016; G. Piano et al. 2017; M. Harvey et al. 2021; Y. Xing & Z. Wang 2021). Previous searches for HE emission from GRS 1915+105 placed only flux upper limits (O. Reimer & A. Iyudin 2003; A. Bodaghee et al. 2013). Intriguingly, excess emission was reported by HEGRA at VHE (F. A. Aharonian & G. Heinzelmann 1998), a hint not confirmed by follow-up observations with either the H.E.S.S. (H. E. S. S. Collaboration et al. 2018) or MAGIC (T. Y. Saito et al. 2009) telescopes.
The relative proximity to Earth of microquasars compared to that of jetted active galactic nuclei (AGNs) makes them a critical tool in the study of particle acceleration in jets. However, as every system detected so far displays a unique phenomenology in the gamma-ray band, the role of microquasars as a population remains ambiguous, especially for LMXBs. The recent detection of two systems in the VHE band suggests that microquasars could provide intermittent contributions to the Galactic cosmic-ray sea at energies up to a few PeV (H.E.S.S. Collaboration et al. 2024; R. Alfaro et al. 2024). However, setting constraints on the magnitude of this contribution remains challenging given the low number of systems detected. Studying the gamma-ray emission from a diverse sample of objects (covering a wide range of jet power, velocity, orientation, and companion mass) is crucial to disentangle the phenomena specific to individual sources from universal properties of particle acceleration in fast-moving jets within our Galaxy.
In this Letter we present evidence for a GeV counterpart to the microquasar GRS 1915+105 using the latest Fermi-LAT data (Figure 1), including almost four times as much data as in the latest study dedicated to the same source (A. Bodaghee et al. 2013). The Letter is structured as follows: Section 2 reports on the observation and analysis details; Section 3 presents the results of the spectral, morphological, and flux variability analyses; Section 4 discusses the possible counterparts and evaluates the feasibility of the association with GRS 1915+105; and Section 5 summarizes our conclusions.
2. Observations and Analysis
We use data from the Large Area Telescope (LAT), a pair-production gamma-ray detector on board the Fermi Gamma-ray Space Telescope (W. B. Atwood et al. 2009). The LAT has been in operation since 2008 August, and thanks to its wide field of view, it surveys the whole sky in photon energies from 30 MeV to more than 100 GeV every 3 hr. The LAT's angular resolution depends strongly on the photon energy, with the point-spread function (PSF) 68% containment radius ranging between several degrees at the lowest energies to less than 01 above 10 GeV (S. Abdollahi et al. 2020). To analyze the LAT data, we use the Fermitools (v2.2.0; Fermi Science Support Development Team 2019) and fermipy (v1.2.0; M. Wood et al. 2017) packages.
We select almost 16 years (from 2008 August 4 to 2024 May 22) of P8R3 SOURCE data (evclass = 128; W. Atwood et al. 2013; P. Bruel et al. 2018) from 1 to 100 GeV within 10∘ of GRS 1915+105. The motivation for the 1 GeV threshold is twofold: to minimize the impact of the bright diffuse emission in our analysis, and to take advantage of the improved PSF at higher energies. Additionally, we apply a cut on the zenith angle at 105∘ to prevent Earth-limb contamination. The selected gamma-ray events are binned into a three-dimensional data cube, consisting of a sky map of dimensions 10∘ × 10∘ (our region of interest; ROI) and 004 bin size, centered at the position of GRS 1915+105 and an energy axis with 16 energy bins equally spaced in logarithmic energy between 1 and 100 GeV.
We fit combined spatial and spectral models to these data using a binned maximum likelihood framework (J. R. Mattox et al. 1996). The analysis is done jointly for the PSF0 to PSF3 event types.4 We use gll_iem_v07.fits and iso_P8R3_SOURCE_V3_PSF3_v1.txt as the Galactic and isotropic diffuse components, respectively, while for the modeling of background sources we employ the 4FGL-DR4 catalog, adding a component for each source (v34; S. Abdollahi et al. 2020; J. Ballet et al. 2023). Energy dispersion correction is applied to all sources except the isotropic diffuse component. In this setup, the detection significance is evaluated considering the test statistic , where is the likelihood value for the null hypothesis and is the likelihood for the tested model. The TS follows a χ2 distribution, and thus TS = 25 approximately corresponds to a detection at 5σ and 4σ for one and four source parameters, respectively.
Given the large source density in the region and the addition of substantial data not accounted for in the 4FGL-DR4 catalog (almost 2 yr), the likelihood fit is performed as follows: (1) we optimize the model of the ROI by fitting the normalization of all sources in descending order of their predicted number of counts (Npred) down to Npred = 1 (through the fermipy.optimize function). Then, (2) we free the normalization of all model components within 5∘ of the center of the region, and all parameters for those sources within 25 or with TS > 500, and fit them with NEWMINUIT. Finally, (3) we search for new sources (>4σ, >02 from known sources) and refit the ROI repeating step 2. These steps are performed twice to stabilize the fit. The procedure results in three new model components close to GRS 1915+105 not present in the catalog, required to explain three (>4σ) gamma-ray excesses in the region. In addition to the TS residual map, we also validate the goodness of our fit through a p-value statistic (PS) data/model deviation estimator (Figures 2(b) and (c); P. Bruel 2021). The distribution of this estimator follows the expected distribution (∝10−a∣x∣, for which we expected a ∼ 1 and find a = 1.08), which allows us to conclude that the region is well modeled. This does not mean that all the new model components not present in the catalog are necessarily of physical origin. For example, they could arise from model defects of the Galactic interstellar emission model (IEM). To address this issue, we perform an additional analysis with gll_iem_v06.fits as an alternative IEM, the diffuse model employed in the latest catalog of hard GeV sources (3FHL; M. Ajello et al. 2017). Of the three excesses, shown in Figure 1 and labeled PS J1915.5+1056, PS J1916.6+1051, and PS J1913.4+1050, only the first is detected with TS > 25 regardless of the IEM, suggesting that PS J1915.5+1056 is the only excess of astrophysical origin.
3. Results
Here we present a study on PS J1915.5+1056, detected with TS = 33.6 with the standard IEM and TS = 47.9 when using the alternative one. It is best described by a pointlike spatial model with best-fit position l, b = (45.409 ± 0.026, −0.292 ± 0.024). Interestingly, it is a potential counterpart for GRS 1915+105, as the known microquasar lies at the edge of the source's containment uncertainty (007 at 99% containment, statistical uncertainty only; Appendix A). Regarding its spectral features, PS J1915.5+1056 is well described with a power law with Γ = 2.25 ± 0.14 (Figure 2(a)) and an integrated photon flux above 1 GeV of (6.74 ± 0.33) × 10−10 photons cm−2 s−1, equivalent to (3.67 ± 0.19) × 10−12 erg cm−2 s−1. We do not find significant evidence for curvature (ΔTS = 2.6) or an exponential cutoff (ΔTS = 3.4) in the spectrum. The quoted uncertainties are statistical only—systematic uncertainties on the flux are mostly dependent on the IEM and of similar amplitude to the statistical uncertainties for standard LAT analyses (F. Acero et al. 2015).
Download figure:
Standard image High-resolution image3.1. Contamination from Nearby Sources
The region surrounding PS J1915.5+1056 is densely populated by gamma-ray sources and possible counterparts other than GRS 1915+105. We discuss associations with sources known via other wavelengths in Section 4.1. Here we discuss instead the possibility that the new excess is due to poor modeling of one of the known nearby GeV sources. In particular, the brightest nearby (<05) source to PS J1915.5+1056 is 4FGL J1916.3+1108, reported to be the GeV counterpart of the supernova remnant (SNR) G045.7-00.4 (H.-M. Zhang et al. 2021). We carry out further morphological tests to determine whether the GeV excess PS J1915.5+1056 can also be attributed to emission from the SNR when accounting for the spatial extension of the GeV emission from the SNR.
We assess the preference among several morphological models (each with k degrees of freedom) by means of the Akaike information criterion (AIC; H. Akaike 1998). This is defined as ΔAIC . Although the improvement is relatively marginal, the best (i.e., smaller ΔAIC) result is obtained by replacing the pointlike spatial model for 4FGL J1916.3+1108 used in the 4FGL-DR4 catalog with a disk, alongside a pointlike component for PS J1915.5+1056 (Table 1). That is, the presence of the new source is preferred even when accounting for the extension of the neighboring source of (ΔTSext = 15.14). This angular size is consistent with the projected size of SNR G045.7-00.4, as well as the results from H.-M. Zhang et al. (2021). If we further test the extension of PS J1915.5+1056, we derive a limit of 014 at 95% confidence (ΔTSext = 4.5). This value increases to 022 if instead we consider the 4FGL J1916.3+1108 as a pointlike source. The preference for an independent source is additionally supported by an analysis restricted to energies between 5 and 500 GeV (Appendix A). Consequently, we can confidently state that the new GeV signal is a distinct source, inconsistent in its spatial and spectral properties with other known nearby sources—notably 4FGL J1916.3+1108 (Γ = 3.00 ± 0.11), even when considering its extension.
Table 1. Likelihood and AIC Values for Different Morphological Analyses
Step | Model | Δk | ΔAIC |
---|---|---|---|
0 | Baseline | ⋯ | ⋯ |
1 | 4FGL J1916.3+1108 (pointlike, 4FGL-DR4 position) | 3 | −129.9 |
2 | 4FGL J1916.3+1108 (pointlike, free position) | 5 | −129.7 |
3 | 4FGL J1916.3+1108 (extended) | 6 | −168.2 |
4 | 4FGL J1916.3+1108 (pointlike) + PS J1915.5+1056 (pointlike) | 9 | −153.0 |
5 | 4FGL J1916.3+1108 (extended) + PS J1915.5+1056 (pointlike, fixed spectrum) | 8 | −170.1 |
6 | 4FGL J1916.3+1108 (extended) + PS J1915.5+1056 (pointlike, free spectrum) | 10 | −169.3 |
Note. The baseline scenario considers the model derived in Section 2 excluding 4FGL J1916.3+1108 and PS J1915.5+1056. The spectrum of 4FGL J1916.3+1108 is always described with a log-parabola shape as in the catalog. The position of 4FGL J1916.3+1108 is fixed to the catalog value in all steps but the second. The position of PS J1915.5+1056 is always fixed to the one obtained in Section 2. We note that, when PS J1915.5+1056 is not included in the model, the best-fit position of 4FGL J1916.3+1108 is always consistent with the reported 4FGL-DR4 position (<009, 99% confidence level) and does not shift in the direction of PS J1915.5+1056 regardless of the spatial model considered.
Download table as: ASCIITypeset image
3.2. Search for Variability
We further explore the temporal properties of the new gamma-ray source, as these could carry additional information about the underlying particle acceleration processes. To this end, we search for periodical signals and flaring episodes, both blindly and specifically to GRS 1915+105. We employ the likelihood formalism developed by M. Kerr (2019) and find no significant periodical modulation, nor any gamma-ray flares (Appendix B). We also find no significant differences between the emission before and after 2018 June, when the system changed its X-ray state. This implies that PS J1915.5+1056 is a persistent gamma-ray source.
4. Discussion
4.1. Identifying the Counterpart
There are many classes of astrophysical objects that could be responsible for a pointlike GeV excess in the Galactic plane. The most frequent association with point sources in the Fermi-LAT data are relatively powerful pulsars (D. A. Smith et al. 2023). We make use of the ATNF Pulsar Catalogue (R. N. Manchester et al. 2005) and search for known pulsars in the vicinity of the excess. There are four pulsars inside a box of 1∘ width around the excess, three of which have too low spin-down power values ( erg s−1) to account for the excess (D. A. Smith et al. 2023). The last object, PSR J1914+1054g, has erg s−1 (S. E. Motta et al. 2023) but lies more than 035 away from the best-fit position (see Figure 3). Given that the derived 99% containment radius for the best-fit position is 007, we can rule PSR J1914+1054g out as a likely counterpart.
Download figure:
Standard image High-resolution imageSNRs can also be bright GeV sources. We query the SNRCat (G. Ferrand & S. Safi-Harb 2012) catalog and find that the only one within the 1∘ box is SNR G045.7-00.4, which lies 03 away from the best-fit position and is already accounted for by existing components in the model. The spectrum of the emission associated with the SNR has a softer spectrum (Γ = 3.00 ± 0.11) than the excess (Γ = 2.25 ± 0.14), which allows us to distinguish them (see Section 3.1). Within an SNR scenario, an asymmetric gamma-ray emission could be accounted for if particles escape from it and interact with nearby dense targets to produce the GeV excess. This has in fact been proposed for another source in the region, 4FGL J1914.5+1107c, whose emission is correlated with dense CO gas at the distance of SNR G045.7-00.4 (H.-M. Zhang et al. 2021). However, the CO map at the distance of the SNR shown in Figure 6 of H.-M. Zhang et al. (2021) reveals no dense target at the location of our GeV excess, ruling out this scenario.
Additionally, several stellar systems within the region could act as particle accelerators. Young (<10 Myr) stellar clusters (M. Cassé & J. A. Paul 1980; M. Ackermann et al. 2011; T. Vieu & B. Reville 2023) are known to emit photons up to the TeV range, although the emission is usually extended. The closest known nearby cluster is Berkeley 43 (E. L. Hunt & S. Reffert 2023), which is not close enough to account for the excess, at more than 03 away from the best-fit position (see Figure 3) with an angular extent of 5'. Besides GRS 1915+105, there is another known X-ray binary system in the ROI, IGR J19149+1036 (G. Cusumano et al. 2015). The GeV excess exhibits no periodicity following the orbital period of IGR J19149+1036 of 22.25 ± 0.05 days (Appendix B). Furthermore, the binary is located more than 035 away from the best-fit position (see Figure 3), so it can be safely discarded as a counterpart.
While the region is filled with possible Galactic counterparts, we also investigate whether the GeV signal could arise from an extragalactic source seen through the Galactic plane. There are no known AGNs in the region in extragalactic catalogs such as NASA/IPAC or WISE (R. J. Assef et al. 2018). However, this is most likely due to the fact that extragalactic catalogs and surveys usually exclude the Galactic plane (D. O. Cook et al. 2023), meaning that the sample of known AGNs at low Galactic longitudes is highly biased to the brightest objects. We used the Fermi-LAT fourth catalog of AGNs (4LAC; M. Ajello et al. 2020) to estimate how many GeV-detectable AGNs are expected in our ROI. There are a total of 2863 objects in the 4LAC catalog, which only considers high Galactic latitudes (∣b∣ > 10∘). This corresponds to a density of 275.7 objects per steradian, leading to 0.08 GeV detectable AGNs expected within the square of 1∘ width shown in Figure 3. Taking into account the spectral properties of the observed excess, the expected number of extragalactic sources goes down by a factor of almost 2 (M. Ajello et al. 2020). We conclude that an extragalactic counterpart is unlikely but cannot at present be fully ruled out.
Finally, the true counterpart could also be an unidentified source. We searched among X-ray catalogs for HE sources other than GRS 1915+105 within the GeV 99% containment region. There are none other than GRS 1915+105 in the SRG/ART-XC all-sky catalog (M. Pavlinsky et al. 2021; S. Sazonov et al. 2024), the INTEGRAL/IBIS 17 yr catalog (R. A. Krivonos et al. 2022), and the second release of the Chandra catalog (I. N. Evans et al. 2024). We note, however, that while the last has a priori the best sensitivity, the imaging coverage of the area was very limited. In the second Swift catalog 2SXPS (P. A. Evans et al. 2020) there are two sources with good detection quality within the excess error radius—both with the same flux and both at the position of GRS 1915+105; thus, they both can be attributed to GRS 1915+105. Finally, in the most recent XMM catalog (N. A. Webb et al. 2020, DR14) there is one nearby source (besides GRS 1915+105) that passes the good-quality flag cut and has a detection likelihood greater than 10. The source is identified as 4XMM J191551.2+105814 (R.A. = 19:15:51.30, decl. = +10:58:14.7), and it has a flux of (1.25 ± 0.19) × 10−13 erg s−1 cm−2 in the 0.2–12 keV energy range, orders of magnitude below that of GRS 1915+105. No variability is reported for this source. The hardness ratio between the 2–4.5 keV and 4.5–12 keV bands is 0.1, indicating a relatively soft spectrum. The much fainter and softer X-ray spectrum than GRS 1915+105 suggests a less powerful source and therefore a less likely counterpart—although, again, this possibility cannot be completely excluded.
We conclude that GRS 1915+105 is the favored counterpart for the new GeV source. In Section 4.2 we discuss a scenario where particles are accelerated in the jets of GRS 1915+105 relatively close to the central binary. Additionally, there have been suggestions that the GRS 1915+105 jets interact with the interstellar medium (ISM) far away from the central binary (L. F. Rodriguez & I. F. Mirabel 1998; C. R. Kaiser et al. 2004; A. J. Tetarenko et al. 2018). There are three regions often invoked in this unconfirmed scenario: the two radio sources (IRAS 19132+1035 and IRAS 19124+1106) and the dense H ii region G45.45+0.06 (all marked by triangles in Figure 3). We do not detect gamma rays at the location of any of these sources.
4.2. Power Requirements
Gamma-ray emission is produced by relativistic particles interacting with their environment. We assume an acceleration region that produces a power-law spectrum of particles dN/dE ∝ E−2 · exp(−E/Ecut) with a cutoff energy Ecut inside the jets of GRS 1915+105 and use the GAMERA code (J. Hahn 2015; J. Hahn et al. 2022) to calculate the GeV emission in both a hadronic and a leptonic scenario. It is unclear whether GRS 1915+105 hosts an alternative source of power such as an equatorial outflow (K. I. I. Koljonen & J. A. Tomsick 2020); therefore, we do not explore this scenario.
If protons are accelerated to relativistic energies in the jets, they would produce GeV radiation owing to inelastic collisions with nearby gas. We use data from the FUGIN 12CO (J = 1–0) survey based on observations at the Nobeyama Radio Observatory (T. Umemoto et al. 2017) to determine the density in the region. The details can be found in Appendix C. We find values ranging between and cm−3 depending on the size of the region considered. For a value of the ambient density of n ∼ 18 cm−3, relativistic protons need to reach energies of at least a few hundred GeV (Ecut ≳ 250 GeV) and have a total energy above 1 GeV of ∼3 × 1049 erg to match the observed gamma-ray flux. If the average power supplied by the jets is 1% of LEdd and at least 10% of that power (a total of ∼1036 erg s−1) goes to acceleration of cosmic rays, supplying that energy would require ∼1 Myr. If the lower value of n ∼ 9 cm−3 is considered, this time goes up by a factor of two, as it is inversely proportional to the density. Given the estimated duration of the current binary evolution stage of GRS 1915+105 of 10.5 Myr (K. Belczynski & T. Bulik 2002), this scenario should in principle be possible, although it is subject to great uncertainty given the unknown real duty cycle of the jets. Once accelerated, the particles would diffuse around the source—how far they can reach depends on the value of the diffusion coefficient D. Assuming, as in D. Khangulyan et al. (2024), an intermediate diffusion regime with D ∼ 1026(E/10 GeV)0.5, 200 GeV protons would have reached scales of around 50 pc over 1 Myr. This corresponds to a 03 radius at the distance of GRS 1915+105. Given the slight offset between PS J1915.5+1056 and GRS 1915+105, it is plausible that particles that escaped from GRS 1915+105 are causing an enhancement in the cosmic-ray density around the source, which then translates to gamma-ray emission only where the target gas is dense enough (F. A. Aharonian 2004; V. Bosch-Ramon et al. 2005; D. Khangulyan et al. 2024).
If instead we assume the emission to be due to electrons, any consideration is subject to the high uncertainty on the magnetic field and radiation fields in the region. Modeling of the radio properties of GRS 1915+105 suggests a high magnetic field on small scales (<1 pc from the engine), on the order of 1 G, which would heavily suppress inverse Compton scattering from electrons (D. Khangulyan et al. 2024). However, the Fermi-LAT results do not tightly constrain the position within the jet. Considering the extension upper limit derived in Section 3, which corresponds to a radius <36 pc at the distance of GRS 1915+105, the gamma-ray emission could arise from a different region where the magnetic field is considerably lower. Assuming a more favorable value of B ∼ 5 μG on the order of the ISM magnetic field (R. Jansson & G. R. Farrar 2012), the interstellar radiation field from C. C. Popescu et al. (2017) at the position of GRS 1915+105, and again an average jet power of 1% of the Eddington luminosity but now with only 1% of this power going to relativistic (>1 MeV) electrons, an injection time several orders of magnitude longer than the estimated age of GRS 1915+105 of 10.5 Myr (K. Belczynski & T. Bulik 2002) would be needed to match the observed gamma-ray flux. One could instead consider an episode in which the accretion rate is higher, which translates to more power carried by the jet and thus available to electrons. Even in such a case, where the jet sustains powers of 10%–100% of LEdd, electrons would require 0.1%–1% of that power for at least 10 Myr to match both the observed flux level and spectral index. The prospects for leptonic emission improve if one considers additional sources of radiation such as the companion star or the accretion disk, with the caveat that these will dominate in regions close to the central engine, where the magnetic field might be much higher as well. Note that the combination of assumed magnetic field of B ∼ 5 μG and interstellar radiation field results in faster loss timescales for inverse Compton than synchrotron for the relevant energies. Higher values of B ∼ 10–20 μG such as those estimated at multiparsec scales from other microquasars (S. Safi-Harb et al. 2022; H. E. S. S. Collaboration et al. 2024) would drive up synchrotron losses and make it even more challenging to explain the GeV observations as inverse Compton emission.
5. Conclusion
We have presented here evidence for a GeV counterpart (PS J1915.5+1056) to the microquasar GRS 1915+105. We find a pointlike source with a relatively hard spectrum, for which we do not find significant transient or periodic flux modulations. The position of PS J1915.5+1056 is inconsistent with other known potential particle accelerators in the region, and its spectral and morphological properties are sufficiently different from those of the GeV counterpart to the nearby SNR G045.7-00.4 to exclude a common origin. Due to a lack of dense gas at the distance of the SNR, we can dismiss the possibility that PS J1915.5+1056 is the result of protons that escaped the SNR. Instead, we find that surrounding GRS 1915+105 there is enough gas to explain the observed emission as arising from the interaction of protons accelerated in the jets with the nearby gas. This scenario relies on a moderate duty cycle of the jet power over a relatively large timescale, which is largely unconstrained and subject to great uncertainty. A leptonic scenario is also plausible but disfavored, as it requires a substantial amount of the jet power, even when assuming a low (5 μG) magnetic field and the interstellar radiation field as targets for interaction. While we cannot fully discard the association with the much fainter X-ray source 4XMM J191551.2+105814 or an unknown blazar seen through the Galactic plane, we note that there is no evidence for a cutoff in the GeV spectrum of PS J1915.5+1056, which suggests that the emission might extend to even higher gamma-ray energies. A detection in the multi-TeV band would rule out an extragalactic origin and further confirm the association with GRS 1915+105, particularly given the improved angular resolution of Cherenkov telescopes. The firm identification of this microquasar as a gamma-ray emitter would establish LMXBs as HE particle accelerators and constrain their contribution to the cosmic-ray content of our Galaxy.
Acknowledgments
We thank Quentin Remy for his help with the gas density estimation, as well as Jian Li, Teddy Cheung, Matthew Kerr, and Philippe Bruel for their comments on the manuscript. This research has made use of the NASA/IPAC Extragalactic Database, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology. This research has made use of the VizieR catalog access tool, CDS, Strasbourg, France (F. Ochsenbein 1996). The original description of the VizieR service was published in F. Ochsenbein et al. (2000). The CO data were retrieved from the JVO portal (http://jvo.nao.ac.jp/portal/) operated by ADC/NAOJ. The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, the K. A. Wallenberg Foundation, the Swedish Research Council, and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase from the following agencies is also gratefully acknowledged: the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Etudes Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515.
Facilities: Fermi-LAT - , No:45m (T. Umemoto et al. 2017).
Software: astropy (Astropy Collaboration et al. 2013), fermipy (M. Wood et al. 2017), godot (M. Kerr 2019), gammapy 1.2 (A. Donath et al. 2023; F. Acero et al. 2024), numpy (Harris et al. 2020), matplotlib (J. D. Hunter 2007).
Appendix A: Analysis Details
The statistical 99% containment region for the best-fit position when using the standard IEM in the analysis above 1 GeV (see Figure 1) is parameterized as an ellipse with semimajor and semiminor axes of 0088 and 0060, respectively, rotated by an angle 143° with respect to the Galactic North (eastward). For the alternative IEM the values are 0095, 0052, and 136°, respectively. A machine-readable version of the flux points shown in Figure 2 can be found in Table 2.
In order to further verify that PS J1915.5+1056 has different properties than the emission from the nearby SNR, we repeat the analysis with a low-energy threshold of 5 GeV (Figure 4). At those energies, PS J1915.5+1056 is still detected at TS = 31, while the TS from 4FGL J1916.3+1108 goes down from TS = 139 to TS = 16. When fitting the region without the additional model component for PS J1915.5+1056, the new best-fit position for 4FGL J1916.3+1108 is found to be at (l,b) = (45.432 ± 0.028, −0.298 ± 0.028) and thus is consistent only with the position of PS J1915.5+1056.
Table 2. Flux Point Values for Both the Standard and Alternative Diffuse Models
Energy | Energy Range | Flux | TS | Flux | TS |
---|---|---|---|---|---|
Standard | Standard | Alternative | Alternative | ||
(GeV) | (GeV) | (10−12 erg cm−2 s−1) | (10−12 erg cm−2 s−1) | ||
1.33 | 1–1.78 | <1.922 | 1.99 | 1.716 ± 0.587 | 8.80 |
2.37 | 1.78–3.16 | <1.642 | 2.63 | 1.022 ± 0.466 | 4.92 |
4.22 | 3.16–5.62 | 0.925 ± 0.401 | 6.17 | 0.830 ± 0.381 | 5.38 |
7.50 | 5.62–10 | 1.289 ± 0.399 | 15.42 | 1.367 ± 0.392 | 18.81 |
13.33 | 10–17.78 | 0.686 ± 0.637 | 5.74 | 0.680 ± 0.344 | 6.41 |
23.73 | 17.78–31.62 | 1.165 ± 0.547 | 9.09 | 1.036 ± 0.505 | 8.61 |
42.17 | 31.62–56.23 | <0.794 | 0 | <1.038 | 0.39 |
74.99 | 56.23–100 | <0.685 | 0 | <0.777 | 0 |
A machine-readable version of the table is available.
Download table as: Machine-readable (MRT)Typeset image
Appendix B: Search for Variability with Photon Weights
We carry out a blind search for state changes by incorporating photon weights into the Bayesian Block algorithm (J. D. Scargle et al. 2013; M. Kerr 2019). Thus, we apply it to the underlying likelihood distribution itself and not to the derived flux points. We select all photons above 100 MeV within 1∘ of the best-fit position of the source, extending the maximum zenith angle to 105∘. The previously derived model ROI is then employed to compute the probability for each photon to be associated with the new source, which is then used as a weight to compute the number of distinct blocks in our light curve (Nb). We penalize the addition of further blocks by means of a prior . Consequently, the parameter γ sets the significance of flux changes. We derive the false-positive rate by redistributing photons randomly among cells (M. Kerr 2019). In order to ensure zero false detections within the 16 yr of data (hence, a false-alarm rate of 10−4), the prior has to be set to γ = 10. From this, we find the data to be well represented by a single cell, and therefore the source does not flare significantly (Figure 5). This is consistent with the results by A. Bodaghee et al. (2013) using 4 yr of data. Additionally, we divide the data set into two ranges: before and after the state change in GRS 1915+105 that occurred on 2018 July 1 (M. Balakrishnan et al. 2021; S. E. Motta et al. 2021). We perform a regular likelihood analysis on each of these ranges, finding no significant distinction between them. Before the state change, the excess has a significance of TS = 23.0 and a flux above 1 GeV of (3.83 ± 0.94) × 10−12 erg cm−2 s−1. After the state change, the values are TS = 11.4 and (3.41 ± 1.19) × 10−12 erg cm−2 s−1 for significance and flux, respectively.
Download figure:
Standard image High-resolution imageDownload figure:
Standard image High-resolution imageWe search for variability that follows the orbital period of GRS 1915+105 by phase-folding the original data set with said period of (P = 33.85 days, T0 = 55458.68 MJD; D. Steeghs et al. 2013). We consider two and four bins of orbital phase and refit the model of the ROI with the standard likelihood analysis above 1 GeV as described in Section 2. We find no flux modulation along the orbit. Finally, we employ the previously derived weights to carry out a blind search for periodical modulations by means of a weighted aperture photometry light curve (Fermi LAT Collaboration et al. 2012). We derive the Lomb−Scargle periodogram, shown in Figure 5, which shows no signal in the power distribution above 3σ. No significant period is found either when using the likelihood formalism from M. Kerr (2019), which considers background periodical fluctuations as well.
Appendix C: Estimation of the Local Density
We use data from the FUGIN 12CO (J = 1–0) survey based on observations at the Nobeyama Radio Observatory (T. Umemoto et al. 2017) to determine the density in the region. We extract the spectra of 12CO in a box with side 04 centered on GRS 1915+105. The resulting spectrum is shown in Figure 6 and shows a peak at VLSR = 26.5 km s−1, which corresponds to a near distance of 1.6 ± 0.4 kpc and a far distance of 10 ± 0.4 kpc (J. Brand & L. Blitz 1993; T. V. Wenger et al. 2018; M. J. Reid et al. 2019; T. Wenger 2021), the latter consistent with the distance to GRS 1915+105 of d = 9.4 ± 0.6 ± 0.8 kpc (M. J. Reid & J. C. A. Miller-Jones 2023). Note that this peak remains even if the spatial extraction region for the spectra is reduced to a box of side 01. We integrate this peak in each pixel of the box, between velocities of 22.7 and 30 km s−1, which results in the map shown in Figure 6. Adopting a CO-to-H2 conversion factor XCO following the functional form in N. Arimoto et al. (1996) and using the same parameters as in CTA Consortium (2024) results in a conversion factor XCO = 8.45 × 1019 s K−1 km−1 cm−2 at the distance of GRS 1915+105. Following the procedure outlined in, e.g., M.-A. Miville-Deschênes et al. (2017), we derive a total mass inside of a region of 01 radius around the best-fit position of the new GeV source of (within the Fermi-LAT PSF at 10 GeV). Assuming that the mass distribution is homogeneous within a sphere, this corresponds to a hydrogen number density of n cm−3. To explore how robust this estimate is, we select a region of the same size around GRS 1915+105, for which the resulting density is n cm−3. If, instead, we consider a larger radius corresponding to the extension upper limit of 022 (Section 3), the resulting density is lower by a factor 2.
Download figure:
Standard image High-resolution image