The prototype X-ray binary GX 339–4: using TeV $\gamma$-rays to assess LMXBs as Galactic cosmic ray accelerators

GX 339–4 has been detected in the optical bands, but the origin of this emission is still not clear. Tetarenko et al. (2020) recently studied its multiwavelength emission and concluded that the optical emission in bright outbursts like the one of 2010 cannot originate exclusively from irradiation of the accretion disc, because unreasonable amounts of energy would be required. Thermal synchrotron emission from the base of the jets could then be considered a good candidate for the optical emission. On the other hand, GX 339–4 shows a flat spectrum in the radio with a spectral break in the IR band that corresponds to the transition of optically thick to optically thin synchrotron emission (Corbel & Fender, 2002; Gandhi et al., 2011). Extrapolating the optically thin IR emission to the X-ray band, significantly underpredicts the optical flux (Maitra et al., 2009; Gandhi et al., 2011; Tetarenko et al., 2019). Hence, if the optical emission originates in the jets, it must come from a different region compared to the IR (Markoff et al., 2003; Corbel et al., 2013).

Reflection features, including a broad iron emission line, are also evident in the X-ray spectrum of GX 339–4 (Nowak et al., 2002; García et al., 2015; Fürst et al., 2015; Parker et al., 2016; García et al., 2019; Dziełak et al., 2019). A jet synchrotron component that is beamed perpendicularly away from the accretion disc is very unlikely to produce significant relativistic reflection (Markoff & Nowak, 2004; Reig & Kylafis, 2021). Furthermore, Uttley et al. (2011) studied the energy-dependent time lags and found that the instabilities in the accretion disc may be responsible for driving the continuum variability on short and longer-than-second timescales. The large time-lags are due to the travel-time between the illuminating region and the disc where the X-rays are reprocessed, and can be only tens of gravitational radii at most. That indicates that the X-ray continuum should be governed by a single component, and a thermal corona close to the black hole could sufficiently explain it (but also see Mahmoud et al., 2019, for a two-component corona).

Based on these results, we approach the modelling assuming the most conservative case for the jet power: that the radio through IR up to the break is self-absorbed synchrotron from the extended jets, the optical emission is synchrotron emission from thermal particles at the base of the jets, and that the X-ray reflecting power-law is from a separate coronal region.

In this work we use archival quasi-simultaneous data to model the multiwavelength spectrum of GX 339–4 from radio to X-rays during the hard state of the 2010 outburst. We use the radio data obtained by the Australian Telescope Compact Array (ATCA) on MJD 55263 (Corbel et al., 2012), IR data obtained by the Wide-field Infrared Survey Explorer (WISE) on MJD 55266 (Gandhi et al., 2011), optical data obtained by the Small & Moderate Aperture Research Telescope System (SMARTS) on MJD 55263, and X-rays from the Neil Gehrels Swift Observatory/X-ray Telescope (Swift/XRT) on MJD 55262 (Corbel et al., 2012) and Rossi X-ray Timing Explorer/Proportional Counter Array (RXTE/PCA) on MJD 55263 (Corbel et al., 2012). We use the 0.5–4.0 keV XRT and the 3–45 keV PCA X-ray data. The IR data are not simultaneous and were obtained 3 days later, but we use them because they show a spectral break crucial for our interpretation of the whole spectrum (see below). There was no significant variability in this time, hence this is a decent assumption to combine these data (Corbel et al., 2012; Corbel et al., 2013; Connors et al., 2019). We also use the upper-limits in the GeV band set by the Fermi/Large Area Telescope (LAT) $\gamma$-ray telescope during the 2010 outburst to further constrain the highest energy regime of the spectrum (Bodaghee et al., 2013). We provide the energy/frequency ranges and the corresponding flux density of all the data we use in Table 1.

We assume that two compact jets are launched by the accreting black hole with jet base radius $R_{0}$. The power injected into the jets in the comoving frame $L_{\rm{jet}}$ defines the number density of the cold (non-relativistic) protons in the plasma at the base of the jets as,

where $\beta_{0,s}\Gamma_{0,s}c$ is the comoving velocity of the plasma in the jet base assumed to be equal to the speed of sound in a relativistic fluid (Falcke & Biermann, 1995; Markoff et al., 2008; Crumley et al., 2017; Lucchini et al., 2022), $\beta={U_{\rm{e}}}/{U_{\rm{B}}}$ is the plasma beta where $U_{\rm e}$ is the energy density of the electrons and $U_{\rm B}$ is the magnetic field energy density.For simplicity, we assume equal number density of electrons and protons, but we discuss the implication of this assumption in Section 5. We further assume that the electron population at the jet base is injected in a thermal Maxwell-Jüttner (MJ) distribution with a peak-energy of $2.23\,k_{\rm{B}}T_{\rm{e}}$.

We vary the plasma beta at the jet base to define the strength of the magnetic field, which scales inversely with distance along the jet $z$. Assuming the electron enthalpy is not significant, we define the magnetisation of the jet as

where $B_{0}$ is the strength of the magnetic field at the jet base. We do not consider any particular magnetic field configuration (toroidal or poloidal) but merely describe the magnetic field by its total strength $B$.

At some distance $z_{\rm{diss}}$ along the jet axis, energy is dissipated into accelerating a fraction of the thermal particles into a non-thermal power-law. We assume that the accelerated particles carry a fixed fraction of the jet power, and in particular, we conservatively fix the power of the non-thermal leptons to be 0.02$L_{\rm jet}$ and of the protons to be 0.05$L_{\rm jet}$.

We allow $z_{\rm{diss}}$ to vary as a fitted parameter and, for the case of the leptonic populations, we assume constant re-acceleration along the jet, but we constrain the proton acceleration to occur only between $z_{\rm diss}$ and $10\,z_{\rm diss}$ in order to limit the required power. Because the most compact part of the jet produces the non-thermal particles, this dissipation region also corresponds to the region where the synchrotron radiation breaks from flat/inverted due to self-absorption, to steep/optically thin. After predictions by Markoff et al. (2001), Corbel & Fender (2002) confirmed that this break typically falls in the NIR band during hard states, and we chose the epoch here because of high-quality observations by Gandhi et al. (2011) that could pinpoint the synchrotron break frequency to be $\rm{4.6^{+3.5}_{-2.0}\times 10^{13}\,Hz}$. To match this frequency, we fix the particle acceleration region at 2600  $r_{\rm{g}}$ from the black hole (see also Connors et al., 2019).

The accelerated particles follow a power-law in energy of the form

In principle, the power-law index $p$ depends on the acceleration mechanism and may differ between electrons and protons, but we choose to use the same for both populations for simplicity.

In equation 3, $E_{{\rm{max}}}$ is the maximum particle energy constrained by energy losses and/or escape. In this work, the maximum electron energy is limited by synchrotron losses and the maximum proton energy is limited by the lateral escape from the jet region. The maximum attainable energy is self-consistently calculated along the jet by equating the characteristic timescales of the losses to the acceleration timescale. The characteristic acceleration timescale $t_{\rm{acc}}=4E/(3f_{\rm{sc}}\rm{ec}B)$ depends on the acceleration efficiency parameter $f_{\rm{sc}}$ that we take to be close to maximum, namely $f_{\rm{sc}}=0.01$ (Jokipii, 1987; Aharonian, 2004). We plot the characteristic timescales versus the particle kinetic energy for the population of the accelerated protons in appendix A.

The fractional number of accelerated particles with respect to the total number of particles $f_{\rm{nth}}$ depends on the acceleration mechanism as well. This number may not be constant along the jet. We parametrize the density of the accelerated particles following:

where $f_{\rm{pl}}>0$ is a free parameter accounting for our ignorance about the exact nature of the dissipation. ($n_{\rm{nth}}$) $n_{\rm{th}}$ is the number density of the (non-)thermal particles. The physical motivation behind such an assumption is the fact that it leads to the characteristic inverted spectrum between radio and optical wavelengths detected in BHXBs (see discussion in Lucchini et al. 2021).

The minimum energy of the accelerated particles depends on the injected distributions in the base. We assume that the minimum energy for accelerated protons is the rest mass energy ($\rm{m_{p}c^{2}}$). This choice is intended purely to limit the number of free parameters; we discuss its implication below. We take the peak of the MJ distribution $2.23k_{\rm{B}}T_{\rm{e}}$ to be the minimum energy of the accelerated electrons. We further define a heating parameter $f_{\rm{heat}}$

The physical motivation behind this assumption is that along with the electron acceleration, some extra heating has been reported by numerical simulations (Sironi & Spitkovsky, 2009; Gedalin et al., 2012; Plotnikov et al., 2013; Sironi et al., 2013; Sironi & Spitkovsky, 2014; Melzani et al., 2014; Crumley et al., 2017). The value of this parameter is not well constrained, but we set it to be $f_{\rm heat}$< 10 (Sironi & Spitkovsky, 2009, 2011; Crumley et al., 2019).

We assume a standard geometrically thin, optically thick accretion disc truncated at some innermost radius $R_{\rm{in}}$ with temperature $T_{\rm{in}}$ (Shakura & Sunyaev, 1973; Frank et al., 2002). We describe the disc luminosity $L_{d}$ in terms of Eddington luminosity $L_{\rm{Edd}}=4\pi{\rm{G}}M_{\rm{bh}}{\rm{m_{p}c}}/\sigma_{\rm{T}}$. We further assume the existence of a hot electron plasma of temperature $T_{\rm{cor}}$, in a spherical region centered on the black hole, normalized by a radius $R_{\rm{cor}}$, and of optical depth $\tau_{\rm{cor}}=n_{\rm e}R_{\rm cor}\sigma_{\rm T}$. These hot electrons upscatter the disc photons to higher energies. We require the existence of such a plasma to be able to model both the X-ray spectrum and properly account for the measured hard timing lags as mentioned in Section 1 (and see e.g., Connors et al., 2019, and discussion below).

Our results with our new lepto-hadronic multi-zone jet model confirm earlier results that stratified jets can self-consistently reproduce the radio-to-X-ray spectrum, together with a thin accretion disc including reflection (Markoff et al., 2001; Zhang et al., 2010; Kylafis & Reig, 2018; Connors et al., 2019; Lucchini et al., 2022). However, compared to earlier works (e.g., Markoff et al., 2005), we can also better reproduce the significantly inverted radio-to-IR spectrum by introducing a decreasing particle acceleration efficiency along the jets (Lucchini et al., 2021). We see however in Table 3 that the parameter $f_{pl}$ controlling this effect cannot be well-constrained by the data, and we can only set an upper-limit.

Apart from particle acceleration, we require significant electron heating of the thermal population(Sironi & Spitkovsky, 2009; Gedalin et al., 2012; Plotnikov et al., 2013; Sironi et al., 2013; Sironi & Spitkovsky, 2014; Melzani et al., 2014; Crumley et al., 2019) to reproduce both the optical and IR bands as jet synchrotron emission. In particular, we find that the scenario where optical emission originates from the jet base and the IR emission originates from the particle acceleration region $z_{\rm diss}$ is consistent with the data. An alternative scenario is that both the IR and the optical emission originate in a hot flow that consists of thermal and non-thermal electrons (Poutanen & Veledina, 2014; Kosenkov et al., 2020), a scenario that better describes the soft states (Kosenkov & Veledina, 2018). Further simultaneous IR-to-optical observations in the hard state would be able to test this scenario, as well as simultaneous polarisation measurements across the entire optical/IR band (although see e.g., Russell & Fender 2008 for measurements prior to the 2010 outburst).

In both the leptonic and lepto-hadronic scenarios the shape of the radio-to-X-ray spectrum of GX 339–4 looks identical and the radiative mechanisms are also the same. The spectral shape is determined primarily by the jet geometry and dynamics, which are similar between the scenarios. However, for the case of the lepto-hadronic models, where we assume equal number density of accelerated electrons and protons, we require much more power injected into the jet base than for the purely leptonic model, which is a well-known issue with hadronic models (see e.g., Pepe et al., 2015; Abeysekara et al., 2018; Kantzas et al., 2020).

To fit the optical emission with thermal synchrotron emission from the base of the jets while the accelerated particles fit the radio-to-IR, we require high electron temperature. This radiation leads to a curved IC spectrum in the soft X-rays, so another component is required to explain the hard power-law. If it can be confirmed that the optical emission is jet synchrotron (via polarisation for instance), then the need for a second component to fit the X-rays will be more robust. For this reason we have added a simple thermal corona model, which together with reflection, can well account for the X-ray spectrum, but is otherwise independent of the jet parameters. In reality, these components should be linked, but it is well known that spectral information alone is often not enough to probe the detailed geometry of the corona, which is the case in our work here as well (see e.g., Del Santo et al., 2008; Droulans et al., 2010; Reig & Kylafis, 2015, 2021; Kylafis & Reig, 2018; Connors et al., 2019; Cao et al., 2021).

When protons are accelerated, the hadronic interactions contribute with additional flux in the $\gamma$-ray regime of the spectrum. For the scenario with a hard proton power-law index of $p_{\rm p}$= 1.7, producing significant TeV flux detectable by CTA would require a non-physical amount of power dissipated into proton acceleration. By constraining the non-thermal proton power to 5% of the jet power, we see that the TeV flux does not exceed the CTA sensitivity (see Fig. 5). A more typical power-law index of $p_{\rm p}$= 2.2 produces even less GeV and TeV flux. In addition, both of these models require strongly matter-dominated outflows even at their launching point ($\sigma\lesssim 0.1$). Such a low magnetisation raises issues of physicality for these models, since the final bulk Lorentz factor of the flow is expected to be on the order of the initial magnetisation $\sigma$ (Komissarov et al., 2007, 2009; Tchekhovskoy et al., 2008, 2009; Chatterjee et al., 2019). Specifically, BHXB jets consistently show at least mildly relativistic velocities of $\Gamma\sim 2-3$ in several systems (Mirabel & Rodriguez, 1994; Fender, 2001; Fender et al., 2004; Casella et al., 2010; Miller-Jones et al., 2012). Such a low initial magnetisation would struggle to explain the bulk acceleration of the flow unless further energy is available by, e.g., thermal pressure. However, numerical simulations show that a jet ’sheath’ forms where the originally Poynting-flux dominated ’spine’ interacts and entrains the surrounding disc wind, resulting in a region with much lower magnetisation (McKinney, 2006; Móscibrodzka et al., 2016; Nakamura et al., 2018; Chatterjee et al., 2019). The instabilities that form along this boundary are expected to be sites of reconnection and particle acceleration (Rieger & Duffy, 2004; Faganello et al., 2010; Rieger, 2019; Sironi et al., 2020). Thus, although our approach is quite simplistic, it would be consistent with the emission occurring along this boundary as suggested by recent radio observations of AGN jets, such as M87 (Hada et al., 2016) or Cen A (Janssen et al., 2021), and GRMHD simulations (e.g. Móscibrodzka & Falcke, 2013; Davelaar et al., 2018). Although BHXB jets cannot be resolved by current facilities, similar scenarios may apply to them since the systems are likely to be governed by the same physical laws (Heinz & Sunyaev, 2003; Merloni et al., 2003; Falcke et al., 2004).

In Fig. 8 and  8 we plot the total distribution of the primary electrons and protons, respectively, integrated along the jets. The MJ-only distribution at the jet base dominates the lower energy regime, with its peak defined by the free parameter $T_{\rm{e}}$ (see Table 3), while the higher energy electrons originate mostly at the first particle acceleration region $z_{\rm diss}$. The shifting of the thermal peak between the two shows the effect of the $f_{\rm{heat}}$ parameter. The fact that the slope is steeper than $p_{\rm e}=2$ indicates that the synchrotron cooling break occurs below $\sim 10^{9}$ eV.

In Fig. 8 we plot the differential number density of the secondary pairs from pp and p$\gamma$ for the lepto-hadronic model with $p_{\rm p}$= 2.2. We also include for comparison the total distribution of the primary pairs of the jets. We note that the secondary pairs from p$\gamma$ are synchrotron cooled and hence their spectrum is flat. The excess of particles around $\sim 10^{12}\,\rm eV$ is responsible for the TeV flux of Fig. 5.

Assuming a maximum value of $f_{\rm{sc}}$= 0.01, we see that the compact jets of GX 339–4 can accelerate CRs up to 100 TeV. Consequently, if this is true and moreover the entire population of BHXBs can accelerate CRs up to 100 TeV, then BHXBs may contribute to the Galactic CR spectrum up to the knee depending on their total number (see also Cooper et al., 2020).

The uncomfortably high proton powers needed for lepto-hadronic jet models has been a topic of discussion for many years (see e.g., Böttcher et al., 2013; Zdziarski & Böttcher, 2015; Liodakis & Petropoulou, 2020; Kantzas et al., 2020). As discussed in Section 5.1, given what we see in AGN jet observations and simulations, we would expect proton acceleration to happen primarily at the interface between the spine and the sheath of the jet, a region of limited volume (Rieger & Duffy, 2019). In our current setup, as a first approximation, we can limit the volume where proton acceleration occurs by reducing the extent of this region with respect to the total jet length. In particular, similar to previous studies (Romero & Vila, 2008; Vila & Romero, 2010; Zhang et al., 2010; Pepe et al., 2015; Hoerbe et al., 2020), we terminate the proton acceleration at a distance 10 $z_{\rm diss}$ from the region where acceleration initiates. As a consequence, we see that even for a hard power law index of $p_{\rm p}$=1.7, the TeV emission of GX 339–4 due to hadronic processes will not be detectable by CTA, but the energy budget remains within reasonable values.

A further way to constrain the total power of the accelerated protons is by increasing the minimum energy of the accelerated particles (Zdziarski & Böttcher, 2015; Pepe et al., 2015). We nevertheless decide to use as the minimum energy for the accelerated leptons the peak of the MJ distribution and for the accelerated protons the rest mass energy (see Section 3), but will explore this in more detailed future work.

Recent high resolution magneto-hydrodynamic simulations have shown that jets can be significantly mass-loaded via instabilities at distances well beyond the launching point (Chatterjee et al., 2019). This progressive mass-loading could significantly reduce the total proton power and make the hadronic models more viable, but this is a project we will pursue in the future.

In both lepto-hadronic models, the optical emission is produced in the jet base due to synchrotron emission from the thermal leptons. The GeV-to-TeV $\gamma$-ray emission on the other hand, is produced in the particle acceleration region and above, which is located at some distance of 3000 $r_{\rm g}$ from the black hole, two orders of magnitude further away from the jet base. Moreover, the $\gamma$-ray is beamed away making it difficult for any attenuation on this optical emission. The IR emission of GX 339–4 is produced in the particle acceleration region where the $\gamma$-ray emission originates as well. We therefore examine any $\gamma$-ray attenuation on the IR emission.

We calculate the optical depth of a 3 TeV $\gamma$-ray that has the maximum likelihood to interact with the $\sim$0.08 eV IR emission using equation 16 of Mastichiadis (2002):

where $\epsilon_{\rm ph}$ is the target photon energy and $n_{\rm ph}$ is the target photon number density of the particle acceleration region. Such a small values indicates that the particle acceleration region is optically thin to TeV $\gamma$-rays.

Despite the fact that GX 339–4 is considered a ’canonical’ low-mass BHXB, it is also amongst the most distant ones. There are Galactic low-mass BHXBs that are as close as approximately 1–3 kpc, e.g., GRO J0422+32 (Webb et al., 2000; Gelino & Harrison, 2003; Hynes, 2005), XTE J1118+480 (Gelino et al., 2006; Hernández et al., 2008), XTE J1650–500 (Homan et al., 2006; Orosz et al., 2004), GRO J1655–40 (Hjellming & Rupen, 1995; Shahbaz et al., 1999; Beer & Podsiadlowski, 2002), GRS 1716–249 (Remillard & McClintock, 2006), GS 2000+251 (Casares et al., 1995; Barret et al., 1996; Harlaftis et al., 1996), V404 Cyg (Miller-Jones et al., 2009), VLA J2130+12 (Kirsten et al., 2014; Tetarenko et al., 2016b), Swift J1357.2–0933 (Shahbaz et al., 2013; Torres et al., 2015), MAXI J1348–630 (Chauhan et al., 2020) and many more at unknown distances that might also be as low as 2–3  kpc (see Liu et al., 2007; Kreidberg et al., 2012; Tetarenko et al., 2016a).

For this reason we also check whether some BHXBs at a distance of 3 kpc with the same $\gamma$-ray luminosity and spectrum as GX 339–4 could be detected by CTA. Assuming that the jets in this putative source have identical properties to GX 339–4, the $\gamma$-ray flux of a nearer source scales as $\left(d_{\rm GX~{}339-4}/d_{\rm source}\right)^{2}\,F_{\gamma}$, where $d_{\rm GX~{}339-4}$ and $d_{\rm source}$ are the distances of GX 339–4 and the source, respectively, and $F_{\gamma}$ is the $\gamma$-ray flux of GX 339–4. We plot this $\gamma$-ray flux in Fig. 9 and compare it to the simulated sensitivity of CTA for various energies, as a function of observation time ¹¹1http://www.cta-observatory.org/science/cta-performance/. The energy range we study here coincides with the energy range of Fermi/LAT which as we can see in Fig. 9 is orders of magnitude less sensitive than CTA for short integration times. We assume that the $\gamma$-ray flux remains constant for up to one day and its uncertainty is of the order of 30 percent. We see CTA is sensitive enough to detect the 100 GeV emission of a GX 339–4-like source at 3 kpc distance, with an exposure of approximately one hour, assuming the emission remains persistent for that long. Consequently, CTA should be able to detect GeV $\gamma$-rays from several future bright outbursts of nearby Galactic BHXBs assuming that the accelerated particles form hard spectra within the relativistic jets produced at peak hard/hard-intermediate states.

We finally examine a more specific example, in particular that of MAXI J1820+070, which is at 2.96 kpc (Gandhi et al., 2019; Atri et al., 2020). During its outburst in 2018, the source was monitored across the multi-wavelength spectrum, from radio to X-rays (Tucker et al., 2018). Here, we merely benchmark the spectral energy distribution instead of optimising to determine the best fit, with the goal of illustrating the similarities and differences with our results on GX 339–4. We use the radio-to-X-ray spectrum, as presented by Tetarenko et al. (2021). We set the black hole mass at 8.5 $\rm{M_{\odot}}$ (Torres et al., 2020), the inclination angle at 63^∘ and the injected jet power at 15% of the Eddington luminosity (Atri et al., 2020). We take the same model parameters we found for the best fit of GX 339–4 for the case of $p_{\rm e}=p_{\rm p}$= 1.7 and present the spectral energy distribution of the 2018 outburst in Fig. 10 and 11. We see that the radio-to-X-ray spectrum is similar to the one of GX 339–4, namely the radio spectrum is due to non-thermal synchrotron radiation, the optical band is due to thermal synchrotron in agreement with Tetarenko et al. (2021) (although see Veledina et al. (2019) for further contributors), and the X-ray spectrum is due to a thermal corona. In contrast to GX 339–4, the p$\gamma$ emission exceeds the CTA sensitivity in the sub-TeV regime. We further compare our predicted spectrum in Fig. 11 to the upper limits set by Fermi/LAT and the Cherenkov telescopes MAGIC, VERITAS and HESS (Hoang et al., 2019). We see that the predicted emission exceeds the upper limits of HESS and marginally those of VERITAS, but it is worth mentioning that these upper limits are derived after 26.9 and 12.2 hours, respectively (Hoang et al., 2019). We are unable to capture the timing signature of the TeV emission with the current version of our model, but we moreover do not know yet whether the high-energy emission of these sources is persistent for up to 20-30 hours (Bodaghee et al., 2013). If MAXI J1820+070’s TeV emission persists for at least a couple of hours during its next outburst, it could then be a possible target-of-opportunity for CTA. Moreover, based on the population-synthesis results of Olejak, A. et al. (2020) and on the recent X-ray observations of Hailey et al. (2018) and Mori et al. (2021), Cooper et al. (2020) estimated that a few thousands BHXBs may reside in the Galactic disc capable of accelerating protons to high energy (also see Fender et al. 2005). If these sources spend approximately 1% of their outburst in the hard to hard-intermediate state (Tetarenko et al., 2016a), then CTA might be able to detect a few tens of BHXBs in its first years of operation.

In our current analysis, we assume equal number density of electrons and protons in the jets, similar to previous studies (Vila & Romero, 2010; Connors et al., 2019). Following this assumption, we derive the jet kinetic power to be $2\times 10^{37}\rm erg\,s^{-1}$. Tetarenko et al. (2021) though suggest that the jets of MAXI J1820+070 cannot be proton dominated and constrain the ratio of protons to positrons to be $\sim 0.6$ otherwise the jet kinetic power, which they estimate to be $6\times 10^{37}\rm erg\,s^{-1}$, may reach 18 times the accretion power. We aim to further study the impact of the pair-to-proton ratio to jet evolution and emission in a forthcoming work.