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

Since the discovery of cosmic rays (CRs) over a century ago, their origin remains an open question. Galactic CRs with energy up to the knee ($10^{15}$ eV) are considered to originate from supernova remnants, but this scenario has recently been questioned due to lack of TeV $\gamma$-ray counterparts in many cases. Extragalactic CRs on the other hand, are thought to be associated with accelerated particles in the relativistic jets launched by supermassive accreting black holes at the center of galaxies. Scaled down versions of such jets have been detected in X-ray binaries hosting a stellar black hole (BHXBs). In this work, we investigate the possibility that the smaller-scale jets in transient outbursts of low-mass BHXBs could be sources of Galactic CRs. To better test this scenario, we model the entire electromagnetic spectrum of such sources focusing on the potential TeV regime, using the `canonical' low-mass BHXB GX 339-4 as a benchmark. Taking into account both the leptonic radiative processes and the $\gamma$-rays produced via neutral pion decay from inelastic hadronic interactions, we predict the GeV and TeV $\gamma$-ray spectrum of GX 339-4 using lower-frequency emission as constraints. Based on this test-case of GX 339-4 we investigate whether other, nearby low-mass BHXBs could be detected by the next-generation very-high-energy $\gamma$-ray facility the Cherenkov Telescope Array, which would establish them as additional and numerous potential sources of CRs in the Galaxy.


INTRODUCTION
Accreting supermassive black holes located at the centres of galaxies are the most powerful engines in the Universe, and some of the most interesting laboratories to investigate the physics of extreme gravity. Of particular importance are those active galactic nuclei (AGN) that exhibit relativistic and collimated jets. The underlying physics that unites the accretion of black holes with the large scale jets is still an unanswered problem. These relativistic jets are considered powerful enough to accelerate particles to very high energy, making them likely a source of extragalactic cosmic rays (CRs) that reach energies of at least 10 19 eV (Hillas 1984;Abbasi et al. 2018;Perrone 2020).
CRs are elementary particles and/or atoms of extraterrestrial origin. The resulting CR spectrum covers ten orders of magnitudes ★ E-mail: d.kantzas@uva.nl in particle energy and shows two very well known characteristic spectral features where the slope changes. The first one is the 'knee' that is located around 10 15 eV, and the second feature is the 'ankle' that is located around 10 17 eV. Current models assume that CRs up to the knee are produced within the Milky Way, while CRs from above the ankle are of extragalactic origin (Hillas 1984;Drury 2012;Blasi 2013). Supernova remnants have long been considered the dominant source of Galactic CRs based on their size and measured magnetic fields (Hillas 1984;Völk et al. 2003;Vink 2012;Ackermann et al. 2013), but due to the lack of TeV -ray counterparts the debate is still open (Aharonian et al. 2019). Given the ability of AGN jets to accelerate cosmic rays, another promising alternative source could be the Galactic jets launched in X-ray binaries comprised of a stellar accreting black hole and a companion star (BHXBs;Mirabel & Rodriguez 1994;Fender 2001;McClintock et al. 2006). Such Galactic jets share the physical properties of AGN jets but on much smaller scales (Heinz & Sunyaev 2002;Romero et al. 2003;Romero & Orellana 2005;Fender et al. 2005;Romero & Vila 2008;Vila & Romero 2010;Vila et al. 2012;Pepe et al. 2015;Cooper et al. 2020;Kantzas et al. 2020).
The presence of jets in low mass BHXBs is transient, and is tightly connected to the properties of the accretion flow. In the so-called hard state, BHXBs display a flat or inverted radio-to-IR spectrum associated with jet synchrotron emission analogous to AGN jets (Blandford & Königl 1979;Hjellming & Johnston 1988;Falcke & Biermann 1995;Markoff et al. 2001;Fender et al. 2006;Corbel et al. 2003Corbel et al. , 2012. BHXBs transit from quiescent to hard and soft states within 'human-like' timescales, hence we can observe the jet launching and jet quenching in real-time (see e.g., Russell et al. 2020). The dynamical timescales are roughly proportional to the mass of the black hole, so it would take typically millions of times longer to detect similar state transitions in AGN.
Accelerated particles in AGN jets are the source of the nonthermal radiation detected over the entire electromagnetic spectrum, from radio to TeV -rays (see e.g. Tavecchio et al. 1998;Celotti et al. 2001;Aharonian 2004;Georganopoulos et al. 2006;Marscher et al. 2008;Ghisellini et al. 2009). However, the exact radiative mechanism has been under debate for a long time because it is tightly connected to the jet composition and the exact particle acceleration mechanism, which remain debated. Two scenarios are generally considered depending on the jet launching mechanism. First, a purely leptonic jet powered by the black hole spin (Blandford & Znajek 1977) may accelerate electrons/positrons that are responsible for the entire multi-wavelength spectrum (Maraschi et al. 1992;Dermer & Schlickeiser 1993;Marcowith et al. 1995;Böttcher & Schlickeiser 1997;Georganopoulos et al. 2002;Ghisellini et al. 2010). Second, a lepto-hadronic jet powered by the accretion disc (Blandford & Payne 1982) (or which starts out leptonic and entrains hadronic mass) may accelerate leptons and baryons that contribute in different energy bands via different mechanisms (Mannheim 1993;Rachen & Biermann 1993;Mücke et al. 2003;Böttcher et al. 2013;Liodakis & Petropoulou 2020).
Recent GeV observations of the high-mass BHXBs Cygnus X-3 (Tavani et al. 2009) and Cygnus X-1 (Tavani et al. 2009;Malyshev et al. 2013;Zanin et al. 2016), and TeV observations of SS 433 (Abeysekara et al. 2018) suggest that some Galactic jets can accelerate particles to high energy. However, it is not known whether all BHXBs, especially the more abundant population of low-mass BHXBs can routinely produce -rays. Until now, only the highmass BHXBs that are characterised by the presence of a strong stellar wind that interacts with the jet have been detected in the GeV and TeV bands (see e.g., Bodaghee et al. 2013).It is thus important to investigate whether the far more populous low-mass BHXBs can also produce -rays. In this paper we approach this question by studying the 'canonical' low-mass BHXB source GX 339-4, extending our previous work on the 'canonical' high-mass BHXB Cygnus X-1 (Kantzas et al. 2020). Similar to AGN jets, the emitting mechanism responsible for any -rays remains unclear, with both leptonic and hadronic processes considered feasible. We are also interested in exploring how the different composition scenarios may affect the jet dynamics and the interpretation of the jet properties.
In this work, we employ a multi-zone jet model to study the hadronic interactions within the jets, as well as the effect on the dynamics and the electromagnetic signature of low-mass BHXB jets. We examine the bright outburst of GX 339-4 in 2010 to model the radio-to-X-ray spectrum with the goal of predicting the TeV radiation originating in the jets. Using the case of GX 339-4 as a model, we assess the likelihood of other, closer low-mass BHXBs to be po-tential sources for the next generation -ray facilities, particularly the Cherenkov Telescope Array (CTA). Such TeV emission may be the signature of efficient CR acceleration inside the BHXB jets, and hence the entire Galactic population of BHXB jets may contribute to the Galactic CR spectrum.
In Section 2 we discuss the physical properties of GX 339-4 and its spectral behaviour. In Section 3 we describe the model we use to study the spectrum of GX 339-4. We present our results in Section 4, discuss their implication in Section 5 and come to our final conclusions in Section 6.

GX 339-4
GX 339-4 is a 'canonical' low-mass BHXB discovered in 1973 (Markert et al. 1973). It undergoes outbursts every two-to-three years that last from a few weeks to months (Belloni et al. 1999;Corbel & Fender 2002;Corbel et al. 2003;Zdziarski et al. 2004;Homan et al. 2005;Belloni et al. 2006;Motta et al. 2009;Corbel et al. 2012). During outbursts, GX 339-4 rises out of quiescence and launches compact jets that contribute to the radio-to-optical spectrum as the source continues into the hard state (Fender 2001;Fender et al. 2004;Corbel et al. 2000Corbel et al. , 2003Corbel et al. , 2012Corbel & Fender 2002;Casella et al. 2010;Homan et al. 2005;Gandhi et al. 2011). Such consistent, repetitive behavior along with extensive and often simultaneous multiwavelength monitoring makes GX 339-4 a perfect target to better understand the properties of relativistic jets.
Although GX 339-4 is a well-studied source, its physical parameters are not well constrained because of the weakness of its companion star. Based on optical photometry the orbital period is estimated to be between 14.8 hours and 16.8 hours (Callanan et al. 1992;Cowley et al. 2002, respectively). The inclination angle is still unknown but is constrained to < 60 degrees because of the lack of eclipsing (Cowley et al. 2002), and the lack of a detection of the companion star means the mass of the black hole is also uncertain. Various current estimates put the mass (in M ) between 4 − 16 (Shidatsu et al. 2011), 5.8 ± 0.8 for an orbital period of 1.75 days (Hynes et al. 2003), > 7 (Muñoz Darias et al. 2008 or 9.8 for a mass function of 1.91 ± 0.08 M (Heida et al. 2017). We adopt the most recent value of bh = 9.8 M of Heida et al. (2017). Hynes et al. (2004) set the distance of GX 339-4 higher than 6 kpc, and Zdziarski et al. (2004) derived a value of 8 kpc while Parker et al. (2016) found a distance of 8 ± 0.9 kpc, which is the distance we adopt here.

Observational constraints in the hard state
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 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.

Observational data
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. , 2013Connors et al. 2019). We also use the upper-limits in the GeV band set by the Fermi/Large Area Telescope (LAT) -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.

MODELLING
In this section we briefly discuss our model, focusing on the interpretation of the free parameters we fit for. A more detailed description of the model can be found in Kantzas et al. (2020) and in Lucchini et al. (2022).

Jet properties
We assume that two compact jets are launched by the accreting black hole with jet base radius 0 . The power injected into the jets in the comoving frame jet defines the number density of the cold (non-relativistic) protons in the plasma at the base of the jets as, where 0, Γ 0, 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), = e / B is the plasma beta where e is the energy density of the electrons and 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 B 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 . Assuming the electron enthalpy is not significant, we define the magnetisation of the jet as where 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 .

Particle acceleration
At some distance 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 jet and of the protons to be 0.05 jet .
We allow 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 diss and 10 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 4.6 +3.5 −2.0 × 10 13 Hz. To match this frequency, we fix the particle acceleration region at 2600 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 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, 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 acc = 4 /(3 sc ecB) depends on the acceleration efficiency parameter sc that we take to be close to maximum, namely 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 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 pl > 0 is a free parameter accounting for our ignorance about the exact nature of the dissipation. ( nth ) 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 (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.23 B e to be the minimum energy of the accelerated electrons. We further define a heating parameter heat e,min = 2.23 heat B e .
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 heat < 10 (Sironi & Spitkovsky 2009Crumley et al. 2019).

Leptonic processes
Following Kantzas et al. (2020), the leptonic radiative processes we take into account are cyclo-synchrotron radiation and inverse Compton scattering (ICS), where the cyclo-synchrotron photons are further upscattered via the synchrotron-self Compton mechanism (SSC) along the jets. Further photon targets for the ICS are the disc photons. We also take into account a precise treatment of pair production due to photon annihilation and pair annihilation to electron-positron pairs (Coppi & Blandford 1990;Böttcher & Schlickeiser 1997).

Hadronic Processes
Accelerated protons interact with the bulk cold protons of the jet and, via proton-proton (pp) interactions, lead to pion production. Neutral pions decay into -rays and charged pions into secondary electrons and neutrinos via the muon decay channel (Mannheim & Schlickeiser 1994). Photomeson interactions between the accelerated protons and target photons (p ) lead to similar distributions of secondary particles. The target photons we consider here are: the thermal radiation of the accretion disc and the non-thermal radiation originating in the compact jet. Finally, we also account for photopair interactions that lead to the formation of pairs, after the inelastic collision between protons and photons. We use the semi-analytical formalism of Kelner et al. (2006) and Kelner & Aharonian (2008) for pp and p interactions, respectively. For the full description of the treatment of the cascades, see Kantzas et al. (2020). For the case of GX 339-4, no photon field is significant enough to attenuate the GeV and TeV emission (also see the discussion below).

Accretion disc and thermal corona
We assume a standard geometrically thin, optically thick accretion disc truncated at some innermost radius in with temperature in (Shakura & Sunyaev 1973;Frank et al. 2002). We describe the disc luminosity in terms of Eddington luminosity Edd = 4 G bh m p c/ T . We further assume the existence of a hot electron plasma of temperature cor , in a spherical region centered on the black hole, normalized by a radius cor , and of optical depth cor = e cor 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).

RESULTS
In this section, we present the results for the best fits of our model to the multiwavelength spectrum of GX 339-4. We explore three different model scenarios: one purely leptonic, and two lepto-hadronic models. For the purely leptonic model, we assume that the nonthermal electrons follow a power-law with = 2.2 (Corbel & Fender 2002;Gandhi et al. 2011). For the two hadronic models, we explore both a soft ( = 2.2) and a hard ( = 1.7) particle power-law, respectively. For all models we fix some common parameters as shown in Table 2. We choose the ratio between the height of the jet base and its radius to be constant and equal to 2 (Maitra et al. 2011;Crumley et al. 2017). The maximum height of the jet is fixed at a large enough value, so it does not influence the spectrum in the radio band via the self-absorption cutoff, and we choose the maximum reasonable particle acceleration efficiency parameter sc = 0.1, which results in maximum proton energies of the order of tens of TeV in the hadronic models. We tie the truncation radius of the thin accretion disc to the jet base radius to reduce model degeneracy because the disc does not contribute to the electromagnetic spectrum at all. We use the tbabs model to account for the neutral photoelectric absorption in the intergalactic medium, using the cross-sections by Verner et al. (1996) and the cosmic abundances by Wilms et al. (2000), where the absorption coefficient H sets the X-ray absorption column. We use the non-relativistic reflect function to treat in a simplified way the reflection detected in GX 339-4, parametrised primarily via the reflection fraction refl = Ω/2 , which indicates the amplitude of the reflected spectrum (Magdziarz & Zdziarski 1995). We choose this simple model in order to minimize the free parameters used to describe the X-ray spectrum, which is well-fit by a power-law. Our focus is on constraining the jet physics that drives the -ray band, thus we retain most of the free parameters for that model.
We R 0 disc innermost radius (see Table 3) out,disc g 10 5 disc outermost radius N H (10 22 cm −2 ) 0.6 absorption coefficient ‡ refl 0.29 reflection fraction Table 2. The fixed parameters of our models, see text for further discussion. Houck & Denicola 2000) to forward fold the model into X-ray detector space, and to find the statistical best fit to the data presented in Section 2.2. We use the emcee function to explore the parameter space using a Markov Chain Monte Carlo (MCMC) method (Foreman-Mackey et al. 2013). We initiate 20 walkers per free parameter and perform 10 4 loops. We reject the first 50% of the run as the "burn-in" period. We provide the 1 uncertainties in Table 3, along with the results of the best fit for each model. In Figures 1-3 we show the best fits of the multi-wavelength spectrum of GX 339-4 for the three different models we explore. In Fig. 1, we show the purely leptonic model, whereas in Figures 2  and 3 the results of the lepto-hadronic models.
The unique contribution of the hadronic processes can only be seen in the TeV -ray band, because the purely leptonic model cannot produce significant emission at GeV and above. In Figures 4  and 5 we show the predicted GeV to 100 TeV -ray spectrum of GX 339-4. The primary-accelerated electrons dominate in the GeV regime via SSC. The hadronic processes dominate in the TeV energy band, in particular, the neutral pion decay from both pp and p collisions as well as the synchrotron radiation of secondary pairs from the latter. Because we set the acceleration efficiency parameter sc to a high value, the protons are able to achieve high energies of the order of ∼ 10 13 eV, producing -rays of the order of TeV. (G) @ diss 1 × 10 4 6 × 10 3 1 × 10 4 p,max (eV) -2.8 × 10 13 2.7 × 10 13 e,max (eV) 1 × 10 8 5.2 × 10 10 5.3 × 10 10 Table 3. Parameters for the three fitted models, distinguished via the power-law index of the accelerated electrons e and protons p . We show the free parameters and the 1 uncertainties as discussed in Section 3 before the double line. Below the double line are indicative evaluated quantities of the plasma magnetisation, the magnetic field, the total luminosity of the accelerated proton/electron population and the maximum energy of the protons/electrons at the particle acceleration region.

Multi-wavelength spectrum and jet dynamics
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 . We see however in Table 3 that the parameter 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 &  sion. 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 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 radioto-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 Xray 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 2015Kylafis & Reig 2018;Connors et al. 2019;Cao et al. 2021).
When protons are accelerated, the hadronic interactions contribute with additional flux in the -ray regime of the spectrum. For the scenario with a hard proton power-law index of 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 = 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 ( 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 (Komissarov et al. 2007(Komissarov et al. , 2009Tchekhovskoy et al. 2008Tchekhovskoy et al. , 2009Chatterjee et al. 2019). Specifically, BHXB jets consistently show at least mildly relativistic velocities of Γ ∼ 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).

Particle distributions
In Fig. 6 and 7 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 e (see Table 3), while the higher energy electrons originate mostly at the first particle acceleration region diss . The shifting of the thermal peak between the two shows the effect of the heat parameter. The fact that the slope is steeper than e = 2 indicates that the synchrotron cooling break occurs below ∼ 10 9 eV.
In Fig. 8 we plot the differential number density of the secondary pairs from pp and p for the lepto-hadronic model with 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 are synchrotron cooled and hence their spectrum is flat. The excess of particles around ∼ 10 12 eV is responsible for the TeV flux of Fig. 4.
Assuming a maximum value of 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).

Non-thermal proton power
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 diss from the region where acceleration initiates. As a consequence, we see that even for a hard power law index of 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.

-ray attenuation on the optical/IR emission
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 -ray emission on the other hand, is produced in the particle acceleration region and above, which is located at some distance of 3000 g from the black hole, two orders of magnitude further away from the jet base. Moreover, the -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 -ray emission originates as well. We therefore examine any -ray attenuation on the IR emission. We calculate the optical depth of a 3 TeV -ray that has the maximum likelihood to interact with the ∼0.08 eV IR emission using equation 16 of Mastichiadis (2002): where ph is the target photon energy and 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 -rays.
For this reason we also check whether some BHXBs at a distance of 3 kpc with the same -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 -ray flux of a nearer source scales as ( GX 339−4 / source ) 2 , where GX 339−4 and source are the distances of GX 339-4 and the source, respectively, and is the -ray flux of GX 339-4. We plot thisray flux in Fig. 9 and compare it to the simulated sensitivity of CTA for various energies, as a function of observation time 1 . 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 -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 -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 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  Figure 9. The -ray light curves in two energy bins as indicated in the legend. The horizontal green lines indicate the predicted flux of a BHXB with the same luminosity as GX 339-4 but located at a distance of 3 kpc instead. We assume that the accelerated particles follow a power-law with index e = 1.7 and p = 1.7, and the emitted flux remains constant for one day.
CTA can detect such a GeV emission within the first hour of the outburst, but Fermi/LAT is not sensitive enough to detect such an outburst.
by (Tetarenko et al. 2021). We set the black hole mass at 8.5 M (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 e = 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 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 highenergy 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 × 10 37 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 ∼ 0.6 otherwise the jet kinetic power, which they estimate to be Flux density (mJy) Syn, z < zdiss Syn, z > zdiss IC, z < zdiss IC, z > zdiss Disc Corona Figure 10. The predicted flux density based on the lepto-hadronic scenario with p =1.7 for the 2018 outburst of MAXI J1820+070. We compare the total emitted spectrum to the data of Tetarenko et al. (2021). The rest of the components are the same as in Fig. 1. 6 × 10 37 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.

SUMMARY AND CONCLUSIONS
Astrophysical jets are ideal laboratories to understand the underlying physics of particle acceleration and the physical processes responsible for the non-thermal emission. It is still unclear whether BHXB jets can accelerate particles to high enough energy to shine in the -ray regime of the electromagnetic spectrum. Such emission strongly depends on the composition of the jets, which remains poorly constrained for either Galactic or extragalactic jets. A possible hadronic composition would support BHXB jets as candidate sources of Galactic CRs and shed light on this long-standing open question. Understanding the jet composition is clearly crucial not only for a better understanding of the non-thermal radiation and total power requirements, but also for our understanding of the jet launching and bulk acceleration properties.
To further understand the properties of Galactic jets and pre-dict any TeV signature, we studied the 'canonical' low-mass BHXB GX 339-4 during the bright outburst of 2010. We presented the best fit of our jet model to the multiwavelength emission and found that the whole radio-to-GeV electromagnetic spectrum can be due to primary leptonic processes. To explain both the radio and the IR/Optical bands, we require a heating mechanism similar to what we see in PIC simulations (Sironi & Spitkovsky 2009Crumley et al. 2019). We further found that the jets of GX 339-4 can accelerate protons to a non-thermal power law up to a few hundreds of TeV. Depending on the power-law index, we saw that the accelerated protons can produce a strong TeV emission via neutral pion decay and synchrotron radiation of secondary pairs. In the case of a hard power law of protons in particular, we found that the photomeson processes dominate the pp interactions and the synchrotron emission of secondary pairs dominates the sub-TeV band. GX 339-4 is however a distant source, located at 8 kpc and the predicted TeV flux will not be strong enough to be detected by future -ray facilities, such as CTA. We rescaled the emitted spectrum to a distance of 3 kpc and compared it to the predicted timing sensitivity of CTA. We find that CTA would be able to detect such emission with an hour of integrated observations in the energy range above 100 GeV, which would be an indication that protons are accelerated into a hard power law. We further tested this scenario by benchmarking the electromagnetic spectrum of a nearby source, such as the newly discovered BHXB MAXI J1820+070. We found that this source might be a potential target-of-opportunity for future CTA observations to hint BHXBs as TeV sources and CR accelerators.