Inverse Compton Emission from Relativistic Jets in Binary Systems

The gamma-ray emission detected from several microquasars can be produced by relativistic electrons emitting through inverse Compton scattering. In particular, the GeV emission detected from Cygnus X-3, and its orbital phase dependence, strongly suggest that the emitting electrons are accelerated in a relativistic jet, and that the optical companion provides the dominant target. Here, we study the effects related to particle transport in the framework of the relativistic jet scenario. We find that even in the most compact binary systems, with parameters similar to Cygnus X-3, particle transport can have a substantial influence on the GeV lightcurve unless the jet is slow, $\beta<0.7$. In more extended binary systems, strong impact of particle transport is nearly unavoidable. Thus, even for a very compact system such as Cygnus X-3, particle transport significantly affects the ability of one-zone models to infer the properties of the gamma-ray production site based on the shape on the GeV lightcurve. We conclude that a detailed study of the gamma-ray spectrum can further constrain the structure and other properties of the gamma-ray emitter in Cygnus X-3, although such a study should account for gamma-gamma attenuation, since it may strongly affect the spectrum above $5\rm\,GeV$.


INTRODUCTION
Microquasars (µQ) are binary systems that host a companion star and an accreting compact object (CO) from which jets are produced. Several microquasars have been detected in the GeV gamma-ray range with AGILE and Fermi LAT (Tavani et al. 2009;Abdo et al. 2009;Williams et al. 2011;Piano et al. 2012;Bulgarelli et al. 2012;Corbel et al. 2012;Malyshev et al. 2013;Bodaghee et al. 2013;Zanin et al. 2016;Piano et al. 2017). The variability found in the GeV emission in some of these sources is consistent with inverse Compton (IC) scattering of stellar photons by relativistic electrons accelerated in the jets (e.g., Dubus et al. 2010a;Zanin et al. 2016;Zdziarski et al. 2018). The IC origin of the gammaray emission detected from µQ is supported by arguments based on higher efficiency of leptonic radiation mechanisms, as compared to hadronic ones, under conditions of compact binary systems (Bosch-Ramon & Khangulyan 2009).
If the dominant target photon field is provided by the stellar companion, IC scattering will be strongly anisotropic E-mail: d.khangulyan@rikkyo.ac.jp (see, e.g., Khangulyan & Aharonian 2005;Khangulyan et al. 2008), and the scattering angle will change along the orbit. This variability of the scattering angle is imprinted in the emission intensity, and may be the dominant factor shaping the GeV lightcurve (e.g., Dubus et al. 2010a, for Cygnus X-3). The specific dependence of the scattering angle on the orbital phase is determined by the jet and counter-jet orientations, and the location of the acceleration and the emission sites in the jet. Thus, gamma-ray light curves can help in constraining the emitter location in µQ.
Cygnus X-3 is the brightest and best studied gammaray emitting µQ (e.g. Tavani et al. 2009;Abdo et al. 2009). The high luminosity of this source may favour, from energetic arguments, relativistic jet velocities, as they could alleviate the demanding energy requirements through Doppler boosting. In such a jet, the non-thermal distribution of particles and their emission would be significantly affected by relativistic effects. Nevertheless, a highly relativistic jet is somewhat in tension with Fermi LAT data in the context of a one-zone IC emitter (Dubus et al. 2010a;Zdziarski et al. 2018). On the other hand, radio VLBI observations of the jets of Cygnus X-3, from milliarcsecond-to-arcsecond scales (∼ 10 − 10 4 AU), favour an at least moderately relativistic jet (Mioduszewski et al. 2001;Martí et al. 2001), which may point to an even more relativistic flow on the scales of the binary (∼ 0.1 AU).
In this paper, we derive the formulas for the IC emission from a relativistic jet using the distribution function of electrons in the phase space (r, p), with r and p being the particle spatial and momentum coordinates in the laboratory reference frame (RF), respectively. This function is a Lorentz invariant, which allows us to avoid cumbersome RF transformations in the case when the contribution from synchrotron self-Compton (SSC) is negligible 1 . This approach also allows us to obtain the results in a form that consistently describes the advection and radiation of gamma rays by particles in the case of an extended emitter. In the derivation, we account both for the transformation of the particle distribution to the laboratory frame, and for the impact of relativistic effects on the particle cooling in the plasma frame. We obtain an analytic solution for the invariant distribution function under the assumption of dominant Thomson IC losses, and numerically compute the IC radiation accounting for changes in the target density and scattering angle along the jet. We discuss the impact of the synchrotron and adiabatic losses, and characterize the conditions when synchrotron losses dominate, under which an analytic solution for the particle distribution can be obtained.
An approach based on the invariant distribution function was earlier suggested to describe the beaming pattern of the external IC emission produced by blobs moving relativistically in blazar jets (Georganopoulos et al. 2001). This approach was later applied to study variable IC emission in binary systems (see, e.g., Kaufman Georganopoulos et al. 2002;Romero et al. 2002). In contrast to these studies, in our paper we consider the emission produced in a jet which implies a different beaming pattern, as compared to an emitting blob. Another difference with the calculations presented by Georganopoulos et al. (2001) is that we use the invariant distribution function to describe the propagation and cooling of relativistic electrons in an extended emitter, which appears to be an important factor for interpreting the gamma-ray emission detected from gamma-ray binary systems.
As compared to other models, which involve extended emitters in µQ (see, e.g., Vila & Romero 2010;Vila et al. 2012;Zdziarski et al. 2014a,b;Pepe et al. 2015) and rely on the conventional approach with RF transformations, our method significantly simplifies the computation of the external IC emission. Thus, this paper allows us to extend the existing models focusing on the GeV gamma-ray emission from µQ (e.g., Dubus et al. 2010a;Zdziarski et al. 2012Zdziarski et al. , 2018, and to study consistently the influence of particle advection on the gamma-ray spectra and lightcurves. Under conditions typical for µQ, in the TeV energy band the Klein-Nishina regime and gamma-gamma attenuation can affect the IC scattering and propagation of gamma rays, respectively (see, e.g., Bosch-Ramon & Khangulyan 2009). In some systems with particularly hot stellar compan-1 Note that in the case of a very clumpy jets, the SSC mechanism may provide a non-negligible contribution (see Zdziarski et al. 2017 ions, e.g. as Cygnus X-3, these effects may influence the production of GeV gamma rays (see, e.g., Protheroe & Stanev 1987;Moskalenko et al. 1993;Bednarek 1997;Cerutti et al. 2011;Sitarek & Bednarek 2012). Therefore, we also consider the influence of the Klein-Nishina effect on the electron transport and the impact of the gamma-gamma absorption on the spectrum adopting system parameters similar to Cygnus X-3.

UPSCATTERING OF EXTERNAL PHOTONS BY A RELATIVISTIC JET
Conventionally the IC radiation from relativistic sources is computed as follows. First one transforms the radiation field to the jet frame. Then, the distribution of high-energy electrons in that frame is obtained, from which the IC emission is computed. Finally, using the relativistic transformation of the radiation to the laboratory frame, one tranforms the emission to that frame (for µQ see, e.g., Dubus et al. 2010b;Zdziarski et al. 2012). Although this method is straightforward, in some contexts it may be more convenient to follow a different path: using Lorentz invariant quantities allows avoiding several frame transformations and provide the results in a form that illustrates the influence of different parameters clearly (see, e.g., Georganopoulos et al. 2001). The difference between these two approaches is sketched in Fig. 1.
To describe the properties of non-thermal particles in the jet, we use the distribution function in the phase space: We follow a general notation policy in which the uppercase "N" refers to number of particles, i.e., a dimensionless quantity, and the lowercase letters, e.g., "n", "f", to densities or (differential) distributions, i.e., dN = ndX, where X is some quantity or a set of quantities. We assume that non-thermal particles are confined in a narrow jet. Thus, in a coordinate system in which the jet is directed along the x-axis, and the system origin is at the CO location, the distribution function should depend on the x-coordinate only, and on the particle momenta. The acceleration site is located at a distance x = x 0 from the CO, and from there non-thermal particles are advected downstream along a relativistic jet that moves with bulk velocity β = V/c and Lorentz factor Γ = 1/ 1 − β 2 (see Fig. 2, where the model geometry is illustrated).
The acceleration site can be associated for instance with a recollimation shock, and may vary with time: x 0 = x 0 (t) (see, e.g., Perucho et al. 2010). However, for the sake of simplicity, it is typically assumed that this distance changes slowly, as compared to the characteristic advection/cooling times, which are also significantly shorter than the binary orbital period (e.g., Bosch-Ramon & Khangulyan 2009).
To compute the radiation accurately, it is necessary to specify the momentum distribution of the particles. The conventional assumption is that the particles are isotropic in the jet co-moving frame, and that are injected following a certain energy distribution (typically a power-law) at x 0 2 . In the jet co-moving frame, the particle distribution function is non-stationary: where x 0 = x 0 /Γ − βct is the location of the acceleration site in the co-moving frame, Θ the Heaviside function, δ the Dirac δ-function, and τ the proper age of the particle population. As non-thermal particles are accelerated at x 0 , one obtains: The evolution of the high-energy particles, which is determined by g , depends also on the parameters that define the jet orientation and x 0 , which for simplicity are assumed to be constant. The internal energy of the plasma per volume unit can be obtained as: 2 The quantities with primes refer to the jet co-moving frame. and the energy injected 3 in non-thermal particles at x = x 0 is: The up-scattering rate of target photons to gamma rays is determined by the differential cross-section: where ε = m 2 e c 4 + p 2 c 2 is the electron energy, with m e being the electron mass, E γ the gamma-ray energy, ε ph the target photon energy, and θ sc the scattering angle, i.e., the angle between the target photon momentum and the observer direction (Aharonian & Atoyan 1981). The scattered photons move in the direction of the electron with accuracy ∼ m e c 2 /ε ( 1 in the relativistic regime), thus an observer located in the direction n 0 should detect emission produced by particles with p = pn 0 .
In the laboratory frame, the gamma-ray spectrum can be obtained from: where θ sc in Eq. (7) depends on r, the location of the source of target photons, and also the orbital phase, the latter playing an important role in shaping the gamma-ray lightcurve. Photons emitted att(r) = T + rn 0 c (T is the time at which a hypothetical photon at r = 0 would be produced), will be simultaneously detected by the observer. The emission will arrive to the observer at a time t = T + T prop , where T prop is the time required for the hypothetical photon to travel to the observer. The linear relation between t and T implies that Eq. (7) describes both emitted and received photons, which is not always the case (see, e.g., Rybicki & Lightman 1979).
The phase-space distribution function is a Lorentz invariant, f (t, r, p) = f (t , r , p ), where t , r , and p are related to t, r, and p by Lorentz transformations (e.g., Landau & Lifshitz 1975). Fort andτ one obtains and where D = (Γ(1 − β cos θ)) −1 is the Doppler boosting factor, with θ being the angle between the jet velocity and the observer direction. The argument of the Heaviside function transforms simply as: Thus, one obtains that the distribution function is stationary in the laboratory frame: but it is anisotropic because of the p x dependence in p = Γ(ε/c − βp x ) − m 2 e c 2 . The integral over dydz can be computed yielding: where is the effective energy distribution density 4 of electrons emitting towards the observer. For the derivation of Eqs. (12) and (13) one considered that the injected particle distribution is relativistic, ε cp Dm e c 2 . Since the target density and scattering angle change along the jet, to obtain the IC emission it is necessary to know the distribution of the particles along the jet.

General case
Since in the jet co-moving frame particles are isotropic, it is convenient to use the energy distribution density: Using the energy-momentum relation, ε 2 = c 2 p 2 +m 2 e c 4 , one can substitute energy as dε = c 2 p dp /ε . Thus, one obtains n = (4π/c 2 )ε p f . For x > x 0 , the energy dependence of the g function can be obtained as For the sake of simplicity, we normalize the particle energy density with Eq. (5), i.e., with a condition that directly determines g . Under the continuous-loss approximation, particle evolution is described as: The energy ε 0 in Eq. (16) evolves accordingly to the particle cooling equation: with the initial condition ε 0 (τ) = ε . The non-thermal energy loss term typically accounts for synchrotron, IC, and adiabatic losses of electrons: in the jet co-moving frame. If the particle distribution function is defined accordingly to Eq. (2), the volume element occupied by the non-thermal particles can be taken as constant dV = dV 0 for a constant-velocity jet (this remains correct independently of the jet geometry -see the discussion after Eq. 29-). Thus, one obtains: If the injection spectrum is a power-law in energy with index α inj , then in the ultrarelativistic limit one obtains: where ε 0 corresponds to ε 0 = ε 0 (0). The normalization constant A is determined by Eq. (5). Then the energy distribution density, Eq. (13), can be represented as for ε = ε/D assuming the relativistic regime, ε ≈ p c.

Compact emitter
When particles are advected along the jet, the intensity of the magnetic and the photon field, and the rate of adiabatic losses, can change, meaning that the above equation does not have an analytic solution in the general case. If advection is slow compared to radiative cooling, Eq. (18) becomes simpler since the loss rates can be considered steady. In this case, the formal solution of Eq. (17) is which should be considered an algebraic equation that determines ε 0 as a function of ε andτ. Once the original electron energy is obtained, one can express the ratio of the infinitesimal energy intervals as

Synchrotron-Thomson losses; extended emitter
The cooling rate for electrons interacting with background photons in the Klein-Nishina regime has a rather complicated dependence on electron energy, which makes finding an analytic solution for Eq. (17) difficult. In contrast, if the IC cooling proceeds in the Thomson regime, the energy loss term has a simple dependence on energy, which is identical to that of synchrotron cooling. One can find an analytic solution if the dominant radiative cooling process is IC in the Thomson regime, or synchrotron emission. In this case, radiation losses have the following dependence on energy: where a does not depend on electron energy. The coefficient a accounts for the synchrotron and IC losses and can vary within the jet, i.e., can be a function of τ: where σ t is the Thomson cross-section; w ph and w b are the energy densities of the target photons and the magnetic field in the jet frame, respectively.
If the dominant photon field is provided by the companion star, then: where L * and R are the luminosity of the star and the distance to it, respectively, and the factor accounts for the transformation of the photon field to the jet co-moving frame (see, e.g., Zdziarski et al. 2012;Khangulyan et al. 2014), where χ is the angle between the jet bulk velocity and the target photon momentum in the laboratory frame (see Fig. 2, where the model geometry is illustrated). The energy density of the magnetic field is: where B is the strength of the magnetic field in the jet comoving frame in G. The rate of adiabatic losses is: where ρ is the plasma density in the flow co-moving frame. We note that ρ corresponds to the jet material. The contribution of the non-thermal particles described by Eq.
(2) to this density might be very small. The shape of the jet (e.g., cylindrical or conical) determines the rate of adiabatic losses. The factor δ(y)δ(z) in Eq. (11) does not imply any limitation on the jet shape. The meaning of this factor is that the conditions for the non-thermal particles do not change considerably across the jet, and Eqs. (2, 11) describe the properties of the particle distribution integrated over the jet cross-section. Equation (18) can be written as: which provides a relation between particle energy at the injection and the emission points as where the relation betweenτ and x comes from Eq. (9). Since the rhs of Eq. (31) does not depend on energy, the ratio of the infinitesimal energy intervals is where ρ 0 and ε 0 are the initial plasma density and particle energy, respectively. The initial energy is For a powerlaw injection spectrum, equation (32) allows us to obtain the electron energy distribution density with Eq. (21). For a non-powerlaw injection, one should use Eq. (19), or equivalently Eq. (15) with the following particle energy distribution: where n 0 is proportional to the injection spectrum.

Synchrotron-Thomson losses; compact emitter
If the rate of radiative losses remains constant over the cooling distance, one can take the parameter a as a constant. However, the plasma density dependence should be preserved in this equation, as otherwise a constant density would imply to ignore adiabatic losses completely. The density evolution is determined by the structure of the jet. For example, for a steady conical jet the mass conservation yields: The evolution of the macroscopic quantities in jets is a subject for dedicated (magneto)hydrodynamic simulations (see, e.g., Perucho et al. 2010), and are beyond the scope of this paper. Equation (35) does not allow determining the density of the jet material, since the jet may undergo bulk acceleration. For the sake of simplicity, we assume that the jet velocity remains constant in the region relevant for gammaray production. In this case, ρ decreases as ∝ x −2 and the integral term in Eq. (34) is: For a short cooling distance as compared to x 0 , this reduces to: which coincides with a solution without adiabatic losses. This calculation illustrates, to some extent, a trivial physical fact: in steady jets the adiabatic losses might be important only in extended emitters. As shown above, assuming a compact emitter (CE) in a steady jet implies a small impact of adiabatic losses. Adopting a power-law in energy with index α inj for the injected particles, one obtains for a compact emitter: Accordingly with Eq. (20), one obtains that in the limit of ε cp Dm e c 2 . The maximum energy in Eq. (39), ε max , is limited by the injection process, and the following relation should be fulfilled: ε /(1 − aτε ) < ε max . In Eq. (12), the limit imposed by the maximum energy translates into the integral upper limit: The dominance of radiation cooling over advection implies that the cooling length, (x max − x 0 ), is small as compared to the characteristic distance over which the loss rate can change significantly. Under these conditions, as noted, one can compute the integral over x analytically: where u max = (ε max − cp/D)/ε max 1, for a large maximum injection energy.

INVERSE COMPTON EMISSION FROM A COMPACT RELATIVISTIC EMITTER
Combining Eqs. (12), (13), (39), and (41), one obtains the following expression for the spectrum produced in a compact gamma-ray emitter in the case of dominant Thomson or synchrotron losses: Accordingly to Eq. (5) for a power-law injection with α inj > 2, the normalization coefficient A is approximately yielding: The relativistic limit, ε ≈ pc, was used above for the sake of simplicity.
If IC scattering of stellar photons is the dominant cooling channel, a ≈ a ic , where a ic is defined by Eqs. (25) and (26), one obtains: where n bb = (2R/R * ) 2 n ph is the Planck distribution, R * the radius of the optical star, and the functionÑ depends on the basic parameters characterizing the system (stellar temperature T * and available non-thermal power), and the acceleration mechanism (accelerated particle energy dependence and minimum energy): where σ b is the Stefan-Boltzmann constant, and the gammaray photon index is related to α inj through α inj = 2(α γ − 1).
If the injection process remains steady over a time similar to, or longer than, the orbital period, there are two factors that affect the variability of the GeV gamma-ray emission. The changing scattering angle and the relativistic effects can vary with the orbital phase. The relativistic effects are accounted by the term D 2α γ +1 D 2 * Γ −1 , which includes both Doppler boosting of the emission, and the transformation of the stellar photon field to the jet RF.
In line with Eq. (13), one can also consider the effective energy distribution of particles in the whole jet for a compact emitter: The emission spectrum can be written as: where the target photon density contains the dilution factor To model the GeV emission from Cygnus X-3, Dubus et al. (2010a) introduced a factor related to relativistic effects obtained for a blob emitter. Zdziarski et al. (2012) argued that this enhancement factor is not applicable to jet sources, and proposed instead the factor derived by Sikora et al. (1997) for the enhancement of the emission in blazars. The impact of stellar field relativistic boosting, D * , is also accounted in the model considered by Zdziarski et al. (2012). In that study, this factor affects several different parameters: electron density; scattering angle; and scattering rate. However, combining their Eqs. (15), (22), and (A9), and taking into account that x (in the notation of Zdziarski et al. 2012) is determined by the scattering angle in the jet frame (x = DD * (1 − cos θ), where θ is the scattering angle in the laboratory frame), one can derive that the IC flux is indeed ∝ D 2 * , which agrees with Eqs. (47), and (48). For a power-law energy distribution of electrons, the gamma-ray spectrum in the Thompson regime allows a simple analytic approximation (Dubus et al. 2010a;Zdziarski et al. 2012). This allows us to obtain the dependence of the IC flux analytically (for the exact expression see Zdziarski et al. 2012):

EVOLUTION OF NON-THERMAL ELECTRONS IN AN EXTENDED EMITTER
Equation (45) describes the IC emission if the cooling length is short as compared to x 0 (see Eq. 40) and the scattering proceeds in the Thomson regime. If ε is the energy responsible for the generation of gamma rays in the laboratory frame, then combining Eqs. (25) and (40) one obtains the condition for a compact production site: with where d is the separation between the normal star and the CO, and fiducial parameter values between those of the gamma-ray emitting microquasars Cygnus X-3 and Cygnus X-1 were used: L * = 10 39 L 39 erg s −1 and d = 10 12 d 12 cm. At typical Fermi LAT energies, say E γ ∼ 200 MeV, the dominant contribution is produced by electron energies of a few GeV (ε =ε GeV). Thus, the condition for the applicability of the compact emitter approximation becomes: and can be violated even for a mildly relativistic jet with Γ ∼ 2.
It is worth noting that for higher energy electrons, IC scattering proceeds in the Klein-Nishina regime. This should lead to an even larger extension of the production site, since IC losses in the Klein-Nishina regime proceed slower than in Thomson.
Losses through IC are not necessarily the dominant cooling mechanism. As indicated above, synchrotron and adiabatic cooling can be also considered. Since adiabatic losses trace a change in the density, they are to be included when the compact emitter approximation fails.
If synchrotron losses dominate, the condition for a compact emitter is determined by the magnetic field in the jet. The synchrotron time is The cooling length, x syn = Γβct syn remains small as compared to x 0 if the magnetic field is: For the sake of simplicity, let us parameterize the jet radius as a fraction of x 0 : r j = x 0 , where 1. Then the Poynting energy flux in such a jet is The numerical coefficient 2.5 × 10 37 corresponds to approximately 5% of the Eddington luminosity for a 5M black hole. Although such a strong magnetic field cannot be excluded in Galactic jet sources, modeling tends to favour a weaker magnetization of the jet (Dubus et al. 2010a;Zdziarski et al. 2012Zdziarski et al. , 2018. Thus, we conclude that there are no robust arguments excluding a significant extension of the gamma-ray production region in gamma-ray emitting microquasars.

Synchrotron and adiabatic loss treatment
The structure of Eq. (31) allows us to consider IC losses and synchrotron losses independently. First we will start with synchrotron losses. The combined impact of synchrotron and adiabatic losses is determined by the following integral: The density and the magnetic field depend on the structure of the jet, but most likely they can be approximated by power-law functions in a limited section of the jet. Thus, the above integral can be solved analytically yielding a description of the joint impact of synchrotron and adiabatic losses. For example, if the impact of IC losses is small, then such a treatment of adiabatic and synchrotron losses allows us to obtain a solution for Eq. (31), and thus to obtain an analytical description of the spatial-energy distribution of electrons in the jet.

Thomson and adiabatic loss treatment
For the case of dominant IC losses, the function a has a more complicated structure since it depends not only on the distance from the CO, but also on the angle between the jet velocity and target photon momenta. Modeling by Dubus et al. (2010a) suggests that the jet in Cygnus X-3 is not perpendicular to the orbital plane. Thus, we consider here also a case when the jet is inclined by a fixed angle α. In this case, the cooling rate depends on the parameter ζ = (x/d−cos α), where α is the angle between the jet velocity and the direction from the CO to the companion star (see Fig. 2, where the model geometry is illustrated). Including ζ, Eq. (17) describing particle cooling becomes: where the dimentionless energy,ε, is defined as In general, Eq. (58) does not allow an analytic solution. We can still treat it numerically, for which we will adopt a density profile of ρ = ρ 0 (x 0 /x) 2 .

Thomson loss treatment
In case adiabatic losses are weak, i.e., ρ is roughly constant, Eq. (58) determinesε(x): where χ = χ(x) and R = R(x) are functions of x, and C is an integration constant. This relation allows one to link the particle energy at the injection point x 0 to that at the emission point x: where Equation (61) together with Eq. (16) allow us to obtain an analytic representation of the energy distribution of the particles in the jet: where the initial particle energy, should remain smaller than the maximum energy, ε max , in the injection spectrum n 0 . For a power-law injection spectrum, one obtains that the distribution function is Equation (65) describes the Thomson cooling of high energy electrons in the jet also in the case when electrons may travel distances comparable to the orbital separation. This analytic solution does not account for the impact of adiabatic losses, as it is the case for instance in a cylindrical jet. Equation (65) can be easily generalized to account for synchrotron losses, for example, if the magnetic field strength has a power-law dependence.
For the computation of the IC emission, one should note that in Eq. (12) the photon density n ph and the IC scattering angle depend on x. In the next section we apply a numerical approach to compute the radiation from such an extended region.

RADIATION FROM AN EXTENDED EMITTER
To study the impact of advection in the jet we adopt parameters similar to those of Cygnus X-3, which is a highly compact system with a very bright star, and also one of the most powerful GeV sources in the Galaxy. The temperature and luminosity of the donor star were adopted to be T * = 10 5 K and L * = 1.8 × 10 39 erg s −1 , respectively, and the CO was assumed to be located at a distance of d = 2.7 × 10 11 cm from the companion star. This separation distance is the largest allowed by Koljonen & Maccarone (2017), and corresponds to a CO with M co 5M and a WR star with M wr 15M .

Emitter size
First, we study the size of the gamma-ray emitting region for different locations of the acceleration site, and for different orientations of the jet. We adopt the maximum energy in the injection spectrum to be ε max = 10 GeV, and compute the maximum distance that electrons with energy ε can reach. For electrons with larger energy, ≥ 10 GeV, and the stellar temperature of Cygnus X-3, the Klein-Nishina effect should weaken the IC losses as compared to the Thomson case. This should lead to an increase of the advection distance, although the impact on the GeV lightcurve should be small given the steepness of the measured spectrum. Efficient advection of high-energy electrons may have an important influence on the multi-GeV gamma-ray flux, since in this energy band the gamma-gamma attenuation is strong, and more efficient advection of the emitting electrons can significantly reduce the attenuation factor. The advection distance, x max , is determined through Eq. (58) as ζ 2 1 + β 2 + sin 2 α − 2βζ ζ 2 + sin 2 α ζ 2 + sin 2 α 2 , which in the case of weak adiabatic losses reduces to 1 We solve numerically Eq. (66) for a conical jet and several jet velocities and orientations. However, if the jet expands slower than a conical one, the rate of adiabatic losses will be smaller. Thus, we also compute the case for dominant IC losses (Eq. 67). Figure 3 shows the relative size of the emitting region, for different injection points: x 0 = 0.3d, d, and 3d (from bottom to top); jet inclinations: α = π/4, π/2, and 3π/4 (from left to right); and jet velocities: β = 0.5, 0.7, and 0.9 (line color). It is seen that even in such a compact system as Cygnus X-3, the advection of GeV electrons might be significant if the acceleration site is located at d ∼ 10 12 cm from the CO, in a jet that moves with a bulk Lorentz factor Γ ≥ 2 (β ≥ 0.87).  . Relative size of the jet filled with electrons of different energies. Calculations are performed for conditions similar to Cygnus X-3, three different jet orientations (columns for α = π/4, π/2, and 3π/4), injection locations (raws: x 0 = 0.3d, d, and 3d), and jet velocities (colors: β = 0.5 -red-, β = 0.7 -green-, and β = 0.9 -blue-). The case with negligible adiabatic losses is shown with dashed lines, and for a conical jet with solid lines.
For a system with a larger star separation than Cygnus X-3, the impact of advection should be stronger. For parameters similar to Cygnus X-1 (T * = 3 × 10 4 K, L * = 8 × 10 38 erg s −1 , and d = 3 × 10 12 cm; see, e.g., Caballero-Nieves et al. 2009, and reference therein), we show in Fig. 4 the extension of the GeV emitter for three different injection points: x 0 = 0.3d, d, and 3d. As seen from the figure, under dominant IC losses the compact emitter approximation for Cygnus X-1 (or similar systems) is only justified if the injection occurs very close to the CO, or if the jet is not relativistic.

Particle distribution
Advection affects the particle energy distribution in the jet. As the baseline case we adopt the results obtained for a compact emitter, Eq. (47). We assume that the nonthermal injection occurs at x = 3d and the acceleration spectrum is steep: α inj = 3. As follows from Eq. (50), for this electron injection spectrum the gamma-ray photon index should be ∼ 2.5, which is roughly consistent with the spectrum detected from Cygnus X-3 with Fermi LAT (see, e.g., Zdziarski et al. 2018). Thomson cooling should render an energy distribution ∝ ε −4 . Thus, in Fig. 5 we plot the particle distributions multiplied by ε 4 . For the case of an extended emitter (EE), we define the effective energy distribution in accordance with Eq. (13) and similarly to Eq. (47) as: In Fig. 5, the results of the calculations of the electron energy distribution are shown for β = 0.5, 0.7, and 0.9, and x 0 = 3d. The maximum energy in the injected spectrum was assumed to be ε max = 30 GeV, the Doppler boosting factor  The case with negligible adiabatic losses is shown with dashed lines, and adiabatic losses for a conical jet with solid lines. D = 1.7 (which implies different viewing angles for different jet velocities), and the magnetic field was assumed to be weak (formally set to B = 0). The energy distribution of a compact emitter, Eq. (47), for the same jet parameters, is shown with gray lines in Fig. 5.
As shown in Fig. 5, the energy distributions that account for advection have three characteristic features as compared to Eq. (47). At the highest energies, there are less particles than in the compact emitter approximation, which is caused by ignoring the injection maximum energy in Eq. (47) (the compact emitter case), and thus it should not be associated with advection. At lower energies, advection leads to particle accumulation in regions with slower energy losses. Thus, for the same injection, the amount of particles for a steady jet is higher. Obviously, these features are the mostly pronounced in electron distributions computed for weak adiabatic losses. Adiabatic losses expected in conical jets appear to be significant enough to determine the shape of the electron distribution at lower energies as shown in Fig. 5 with thin lines.
To illustrate the effect of photon target dilution, we also compute the amount of emitting electrons weighted by the target photon density (thin lines in Fig. 5): As seen in the figure, the weighted energy distribution is suppressed with respect to the compact emitter approximation at low energies because of particle escape. In the calculations, we adopted a maximum jet length of 35d, as particles reaching that far from the CO are already strongly cooled down due to adiabatic losses. We also note that particles that reached that distance from the CO should not produce variable gamma-ray emission, as θ sc does not significantly change along the orbit. Figure 5 shows that for β ≥ 0.7, the particle energy distribution features a break or   (70)). Other model parameters were set as B = 0, x 0 = 3d, and α = π/2. steepening at energies around ε 1 GeV, which might be testable with Fermi LAT in Cygnus X-3.

Spectral energy distribution and gamma-ray lightcurve
To compute the radiative output from a jet in a binary system, it is necessary to define the inclination of the orbit and the jet orientation. The jet orientation should also include the assumption if one considers jet or counter-jet; the former is located on the same side as the observer, and the latter on the opposite side, with respect to the orbital plane. We consider two cases: i orb = 30 • and i orb = 60 • , and the emission produced only by the jet since for relativistic bulk velocities the emission from the counter-jet is strongly suppressed. To illustrate the impact of advection we consider the simplest case of a circular orbit, which is reasonable for Cygnus X-3. For simplicity, the jet is assumed to be perpendicular to the orbital plane. In this case, there are three relevant orbital phases: superior/inferior conjunction of the CO (SUPC/INFC); and when the CO crosses the plane of the sky (NODE). Inverse Compton losses and radiation are affected by three angles: χ (jet velocity-photon momentum angle); θ (jet velocity-line of sight angle); and θ sc (IC scattering angle), which are shown in Fig. 2. These angles depend on the inclination of the orbit, the jet orientation, and the emitter location in the jet. The dependence of these angles on the distance from the CO is shown in Fig. 6 for the selected inclinations. The combined effects of particle advection and changes in the photon density and scattering angle may affect the gamma-ray spectrum in a quite complex manner, as shown in Fig. 7. For certain orientations of the jet, the emission from the inferior conjunction is strongly suppressed by unfavorable scattering angles, in which case advection along x/d θ sc (supc) θ sc (infc) θ sc (node) Figure 6. Characteristic angles that determine electron cooling and emission in a relativistic jet. Figure shows χ (target-photon momentum -bulk velocity), θ (observer direction -bulk velocity), and θ sc (observer direction -target-photon momentum) angles. Calculations are performed for two orbital inclinations: i orb = 60 • (top panel) and i orb = 30 • (bottom panel). The jet was assumed to be perpendicular to the orbital plane, α = π/2, and thus only θ sc depends on the orbital phase; three different orbital phases are shown: SUPC (solid lines), INFC (long dashed lines), and NODE (dashed lines). the jet will tend to enhance the emission below 1 GeV. At high electron energies, the cooling length is shorter, and the IC cross-section has a weaker dependence on the scattering angle, so the particle transport has a minor impact on the high-energy part of the spectrum. Thus, advection may distort the spectral the spectral shape in the GeV energy band. This effect is strongly pronounced in fast jets (see top and bottom panels of Fig. 7).
Around the superior conjunction of the CO, advection tends to harden the radiation spectrum. This is caused by the escape of lower energy electrons to regions of weaker photon fields, as illustrated in Fig. 5 with thin lines. As shown in Fig. 7, advection tends to enhance the emission around inferior conjunction, and to reduce the emission around superior conjunction. Thus, the orbital variation of the gammaray emission will be weakened by advection. This is illustrated in Fig. 8, where lightcurves for different jet velocities (β = 0.7 and β = 0.9), and inclinations (i orb = 30 • and i orb = 60 • ), are shown, with the acceleration site assumed to be at x 0 = 3d. It can be seen that even in the slower case with β = 0.7, the orbital variation becomes significantly weaker.
As expected, adiabatic losses result in a less extended production site, and thus the orbital phase dependence remains stronger than in the case of negligible adiabatic losses (solid vs dashed lines in Fig. 8). However, the spectral change caused by propagation effects may remain strong. In particular, the combined effects of particle advection towards regions with weaker photon field, and adiabatic cooling, result in a hardening of the photon spectrum around 0.1 -1 GeV (see Fig. 7). . Spectral energy distributions of the IC emission for an extended emitter obtained for three different orbital phases: SUPC (green); INFC (red); and NODE (blue). The jet is assumed to be perpendicular to the orbital plane, α = π/2. The case with weak adiabatic losses is shown with dashed lines, and adiabatic losses for a conical jet with solid lines. The region between these two regimes is filled with color. Gray lines show the spectral energy distributions obtained for a compact emitter (Eq. 45). Gray dashed lines show the spectral energy distributions for SUPC and NODE phases obtained by applying the orbital-phase dependent coefficient in Eq. (50) to INFC spectrum obtained under CE approximation (gray solid line). Cases with β = 0.9 and i orb = 30 • , β = 0.7 and i orb = 30 • , and β = 0.9 and i orb = 60 • , are shown in the top, middle and bottom panels, respectively. Other model parameters were set as B = 0 and x 0 = 3d. The jet was assumed to be perpendicular to the orbital plane, α = π/2. The case with weak adiabatic losses is shown with dashed lines, and adiabatic losses for a conical jet with solid lines. The region between these two regimes is filled with color. Two different jet velocities are shown with different colors: β = 0.7 (blue); and β = 0.9 (red). Green lines correspond to the electron spectrum obtained for a compact emitter(Eq. 47). Other model parameters were set as B = 0 and x 0 = 3d.

QUANTUM ELECTRODYNAMICS EFFECTS
There are two quantum electrodynamics (QED) effects that may have a substantial influence on the gamma-ray emission produced in compact binary systems. The first is related to the transition from the classical Thomson limit to the quantum Klein-Nishina regime. The Thomson limit is valid when electron and target photon energies are small: If the electron and the target photon energy are high enough to violate this relation, the precise QED cross-section should be used (for astrophysical conditions, see Aharonian & Atoyan 1981). The Klein-Nishina effect has a strong impact both on the energy loss rate and on the IC spectrum. The second important QED effect is the gamma-gamma attenuation. Typically, in binary systems this effect is important in the TeV energy band (Dubus 2006). If the stellar temperature is high, as, e.g., in Cygnus X-3, the attenuation might be important for gamma rays with relatively low energy, E γ ≥ 10 GeV (Protheroe & Stanev 1987;Moskalenko et al. 1993;Bednarek 1997;Cerutti et al. 2011;Sitarek & Bednarek 2012). In Fig. 9 we show the attenuation factor for gamma rays interacting with the stellar field. The target photon field is provided by the optical star with radius and temperature of R * = 1.6 × 10 11 cm and 10 5 K, respectively. The calculation takes into account the finite size of the star integrating over the stellar surface (which can give a substantial difference as compared to calculations adopting a point-like approximation for the star if the gamma-ray emitter locates within a few stellar radius distance from the star, see, e.g., Dubus 2006;Bosch-Ramon et al. 2008;Bosch-Ramon & Khangulyan 2009;Romero et al. 2010); this accounts for the occultation by the star, as seen in the map opacity for 1 GeV photons in Fig. 9. For 10 GeV gamma rays the attenuation can be very significant for almost half of the orbit unless the gamma-ray production site is located at a large distance from the CO. At a few GeV the influence of the gamma-ray absorption is smaller, although it still can suppress the emission from the counter-jet at SUPC phases, which can be relevant if the jet bulk velocity is relatively small. This may result in an additional factor affecting the orbital variability with a strong energy dependence, yielding a multi-GeV lightcurve significantly different from the GeV lightcurve.
In the case of Klein-Nishina losses, electrons may lose a significant fraction of energy in a single interaction, which is inconsistent with the assumptions used for the continuousloss approximation. However, the continuous-loss approximation was shown to provide results consistent with a detailed kinetic treatment (see, e.g., Khangulyan & Aharonian 2005). Thus, to account for the Klein-Nishina effect we solve Eq. (18) for a conical jet using the approximation for IC losses in a Planckian photon fiedd suggested in Khangulyan et al. (2014). To compute the electron density, we use the continuous-loss approximation, i.e., the energy distribution density of electrons is described by Eq. (21). The influence of the reduction of energy losses due to the Klein-Nishina effect is shown in Fig. 10. It is seen that the weakening of the IC energy losses results in a ∼ 30% increase of the number density of GeV electrons.
For each considered gamma-ray energy and location in the jet we compute the gamma-gamma opacity, τ γγ , in the stellar photon field in the direction of the observer. The influence of different transport and cooling assumptions is illustrated in Fig.11. As seen in the figure, Klein-Nishina IC cooling has a similar impact at different orbital phases on the gamma-ray emission intensity, resulting in a small transformation of the lightcurve shape. In contrast, gamma-gamma attenuation strongly affects the gamma-ray spectrum above 5 GeV for SUPC. We note that Fig. 11 shows the emission produced in the jet; for the counter-jet the impact should be considerably stronger, for relatively low jet velocities. Thus, a detailed study of multi-GeV gamma-ray emission from Cygnus X-3 may significantly constrain the possible locations of the production site. Figure 9. Gamma-gamma attenuation factor for 1 GeV (top panel), 5 GeV (middle panel), and 10 GeV (bottom panel) gamma rays traveling towards an observer looking from the right. The calculations are done for R * = 1.6 × 10 11 cm and T * = 10 5 K.

CONCLUSIONS
We have studied the properties of the gamma-ray emitting region in a relativisitc jet in a binary system. To facilitate the interpretation of the results, we have used an approach based on the distribution function in the phase space, which is Lorentz invariant. This allows obtaining results in a compact form that permits studying the influence of different . Energy distribution of electrons in the jet calculated for conditions similar to those in Cygnus X-3. The jet velocity was assumed to be β = 0.9 and nonthermal electrons with α inj = 3 were injected at a distance x = 3d from the CO. Gray lines correspond to densities shown in Fig. 4 and obtained for Thomson losses: compact emitter (CE, thin line); extended emitter with Thomson losses only (dashed line); and extended emitter with adiabatic losses for a conical jet (solid line). The solid blue line shows the energy distribution of electrons computed using an accurate IC loss prescription, under adiabatic losses in a conical jet. The jet was taken perpendicular to the orbital plane, α = π/2. parameters in a clearer way. The main focus of the study was on the impact of advection on the gamma-ray spectrum and lightcurve.
For the case of a compact production site we have obtained an analytic representation of the energy distribution of the emitting electrons. When IC cooling dominates over advection, the gamma-ray spectrum, given by Eq. (45), has a simple form that allows one to determine the process that affects the variability of the emission. Namely, it contains three factors that change with orbital phase: (i) IC proceeds in the anisotropic regime, and the scattering angle varies along the orbit Dubus et al. 2010a;Zdziarski et al. 2012); (ii) the Doppler boosting factor, D 2α γ +1 Γ −1 , which accounts for the relativistic transformation of radiation produced in a stationary jet (Sikora et al. 1997);and (iii) in the case of dominant IC losses, an additional factor, D 2 * , should be introduced. The stellar photon boosting effect on cooling can be ignored if the dominant losses are due to synchrotron cooling.
Adiabatic losses can be relevant only if relativistic particles are advected along the jet over a distance in which the jet material density undergoes a significant change. In particular, this can be the case for low-energy electrons that are subject to slower radiative losses. In the case of a (at least) mildly relativistic jet, Γ ≥ 2, advection might be important for GeV emitting electrons even in the most compact binaries like Cygnus X-3. In the case of dominant radiative losses, we have obtained an analytic solution that describes the properties of non-thermal electrons in a relativistic inclined jet. This solution can however be generalized to the case when adiabatic losses are important under weak IC losses, i.e., covering a broad range of synchrotron and adiabatic losses. β=0.9, i orb =π/3, x=0.1d Figure 11. Comparison of gamma-ray spectral energy distributions calculated under different assumptions on transport and cooling for: emission obtained for a compact emitter (gray lines); dominant Thomson cooling (green lines); Thomson cooling with adiabatic losses in a conical jet (blue); and Klein-Nishina cooling with adiabatic losses in a conical jet (red lines). Two orbital phases are shown: SUPC (solid lines) and INFC (dashed lines). Calculations account for gamma-gamma attenuation in the photon field of the optical companion. For SUPC, the intrinsic spectra are shown with thin lines (for INFC the attenuation is negligible). The jet velocity was assumed to be β = 0.9, and two different injection points are considered: x = 3d (i orb = π/6, top panel) and x = 0.1d ( i orb = π/3, bottom panel). The jet was taken perpendicular to the orbital plane, α = π/2. It is generally expected that in gamma-ray emitting µQ IC losses should dominate over synchrotron for GeV electrons (see, e.g., Zdziarski et al. 2012). Thus, as test cases, we have considered two cases for extended emitters: (i) dominant IC losses, which allow an analytic solution for the particle density, Eq. (63); and (ii) the case with IC and adiabatic losses, the latter being expected in a conical jet, for which a numerical treatment has been applied. The simulations have shown that, in systems similar to Cygnus X-3, particle advection may have a significant impact on the gamma-ray lightcurve if the jet velocity is high, β ≥ 0.7. For even faster jet velocities, β ∼ 0.9, one should also expect a strong transformation of the gamma-ray spectrum from different orbital phases. In a more extended system, e.g., in Cygnus X-1, advection is very important unless synchrotron losses prevent efficient particle transport (see Eq. (55)), which is probably not very realistic.
In the specific case of Cygnus X-3, the stellar companion should be very hot, T * 10 5 K. For such a target photon field, two QED effects may influence the electron transport and gamma-ray spectrum in the GeV energy band. The Klein-Nishana effect weakens the IC energy losses and affect the gamma-ray spectrum, and gamma-gamma absorption can significantly suppress the flux above a few GeV. To study the influence of these effects we have performed detailed calculations of the electron transport, radiation, and gammagamma opacity. Since in the case of Cygnus X-3, the orbital separation is comparable to the stellar radius, in the calculations of the gamma-gamma opacity we have accounted for the finite size of the optical star. The simulations show that the Klein-Nishina effect has a small impact on the intrinsic gamma-ray spectra. Unless the gamma-ray production site is located at large distance from the CO, x 3d, the gammagamma attenuation should significantly affect the spectrum at multi-GeV energies, E γ > 5 GeV.
To summarize, we have performed a detailed study of the IC process in realistic jets in compact binary systems. The performed study has revealed that the particle advection along the jet might be important even in a very compact binary system, e.g., in Cygnus X-3. In systems similar to Cygnus X-1, advection should be accounted for even in the case of a weakly relativistic jet.
If adiabatic losses are weak, which would be the case, e.g., in cylindrical jets, advection can impact significantly the gamma-ray emission, potentially leading to a strong dependence of the gamma-ray spectrum shape on the orbital phase. For advection in a conical jet, adiabatic losses weaken the effects on the spectrum.
Independently of the dominant cooling channel, advection results in a significant weakening of the orbital phase dependence. Thus, if the properties of the accelerator in Cygnus X-3 and Cygnus X-1 are similar, one should expect differences in the orbital phase dependency of the GeV emission between these two systems.
To illustrate the relevance of this effect, in Fig. 12 we show the lightcurves computed for a system similar to Cygnus X-1 (the temperature and luminosity of the optical star are taken as T * = 3 × 10 4 K and L * = 8 × 10 38 erg s −1 , respectively; the CO was assumed to be in a circular orbit with d = 3.2 × 10 12 cm).
The IC emission shown in Fig. 12 was averaged over two orbital phase bins: |φ| < 0.25 and |φ| > 0.25, the orbit being −0.25 < φ < 0.75). The injection point was assumed to be located at x 0 = 4d, and the injection spectrum and jet velocity were assumed to be ∝ ε −4 and β = 0.5, respectively. The orbital inclination was selected to be i orb = π/3 and the figure includes only the contribution from the jet (i.e., the counter-jet emission is not accounted for because it is expected to be relatively small and to weaken the orbital phase dependence even stronger). The data points are from Zdziarski et al. (2017), and the open and filled squares correspond to the emission expected from a CE and a EE, respectively. The adiabatic losses were assumed to be weak and the magnetic field set to B = 0, so the dominant cooling mechanism is the Thomson scattering. As seen from Fig. 12, the advection may provide a possible explanation for a weaker orbital phase dependence of the GeV emission from Cygnus X-1, and alleviate the requirement for a SSC contribution. Zdziarski et al. (2017) studied the broadband emission and gamma-ray variability in Cygnus X-1 for the parameter space with x 0 d. In that parameter space, it was found that external Compton models, including those with an extended emitter, are incompatible with the Fermi LAT emission from Cygnus X-1, and a better agreement can be achieved if one assumes a highly clumpy jet, which enhances the SSC emission.
In the case of a conical jet (as was assumed by Zdziarski et al. 2017), adiabatic losses lead to considerable cooling at distances x ∼ x 0 , so plasma cools down on a scale in which the IC regime does not change for x 0 d. Thus, advection cannot considerably affect the IC lightcurve for parameters adopted by Zdziarski et al. (2017). The simulations shown in Fig. 12 show that for x 0 ≥ d, advection can improve the agreement between the observational data from Cygnus X-1 and predictions of models that account for stellar IC only. We note, however, that the calculations presented in Fig. 12 are for illustrative purposes only and cannot substitute a detailed broadband study (as the one presented in Zdziarski et al. 2017).
The influence of advection on the gamma-ray light curve also significantly affects the ability of one-zone models (Dubus et al. 2010a;Zdziarski et al. 2012Zdziarski et al. , 2018 to accurately infer the properties of the gamma-ray production sites even in the case of the most compact binary systems. For example, Zdziarski et al. (2018) suggested that the Fermi LAT emission in Cygnus X-3 is best explained by IC scattering from a production site located at x = 2.3d in a jet with β = 0.73. As shown by our simulations, for this location of the production site the transport effects might be relevant.
We present in this paper the theoretical framework and discuss the impact of advection on the GeV gamma-ray spectrum and lightcurve. A detailed application to the gammaray data of Cygnus X-3 will be presented in a forthcoming paper. Cyg X-1 Figure 12. Lightcurves of the IC emission from a system similar to Cygnus X-1: T * = 3 × 10 4 K, L * = 8 × 10 38 erg s −1 ,d = 3.2 × 10 12 cm (circular orbit), β = 0.5, and i orb = 60 • . The jet was assumed to be perpendicular to the orbital plane, α = π/2. The case with weak adiabatic losses is shown with filled squares and the emission expected from a CE is shown with open squares. The data points are adopted from Zdziarski et al. (2017). Other model parameters were set as B = 0 and x 0 = 4d.