Pulsar wind zone processes in LS 5039
Abstract
Several ray binaries have been recently detected by the HighEnergy Stereoscopy Array (H.E.S.S.) and the Major Atmospheric Imaging Cerenkov (MAGIC) telescope. In at least two cases, their nature is unknown. In this paper we aim to provide the details of a theoretical model of close ray binaries containing a young energetic pulsar as compact object, earlier presented in recent Letters. This model includes a detailed account of the system geometry, the angular dependence of processes such as KleinNishina inverse Compton and absorption in the anisotropic radiation field of the massive star, and a Monte Carlo simulation of leptonic cascading. We present and derive the used formulae and give all details about their numerical implementation, particularly, on the computation of cascades. In this model, emphasis is put in the processes occurring in the pulsar wind zone of the binary, since, as we show, opacities in this region can be already important for close systems. We provide a detailed study on all relevant opacities and geometrical dependencies along the orbit of binaries, exemplifying with the case of LS 5039. This is used to understand the formation of the very highenergy lightcurve and phase dependent spectrum. For the particular case of LS 5039, we uncover an interesting behavior of the magnitude representing the shock position in the direction to the observer along the orbit, and analyze its impact in the predictions. We show that in the case of LS 5039, the H.E.S.S. phenomenology is matched by the presented model, and explore the reasons why this happens while discussing future ways of testing the model.
keywords:
rays: theory, Xray binaries (individual LS 5039), rays: observations1 Introduction
Very recently, a few massive binaries have been identified as variable veryhighenergy (VHE) ray sources. They are PSR B125963 (Aharonian et al. 2005a), LS 5039 (Aharonian et al. 2005b, 2006), LS I +61 303 (Albert et al. 2006, 2008a,b), and Cyg X1 (Albert et al. 2007). The nature of only two of these binaries is considered known: PSR B 125963 is formed with a pulsar whereas Cyg X1 is formed with a black hole compact object. The nature of the two remaining systems is under discussion. The highenergy phenomenology of Cyg X1 is different from that of the others. It has been detected just once in a flare state for which a duty cycle is yet unknown. The three other sources, instead, present a behavior that is fully correlated with the orbital period. The latter varies from about 4 days in the case of LS 5039 to several years in the case of PSR B125963: this span of orbital periodicities introduces its own complications in analyzing the similarities among the three systems.
LS I +61 303 shares with LS 5039 the quality of being the only two known microquasars/ray binaries that are spatially coincident with sources above 100 MeV listed in the Third Energetic GammaRay Experiment (EGRET) catalog (Hartman et al. 1999). These sources both show low Xray emission and variability, and no signs of emission lines or disk accretion. For LS I +61 303, extended, apparently precessing, radio emitting structures at angular extensions of 0.010.05 arcsec have been reported by Massi et al. (2001, 2004); this discovery has earlier supported its microquasar interpretation. But the uncertainty as to what kind of compact object, a black hole or a neutron star, is part of the system (e.g., Casares et al. 2005a), seems settled for many after the results presented by Dhawan et al. (2006). These authors have presented observations from a July 2006 VLBI campaign in which rapid changes are seen in the orientation of what seems to be a cometary tail at periastron. This tail is consistent with it being the result of a pulsar wind. Indeed, no large features or highvelocity flows were noted on any of the observing days, which implies at least its nonpermanent nature. The changes within 3 hours were found to be insignificant, so the velocity can not be much over 0.05. Still, discussion is ongoing (e.g. see Romero et al. 2007, Zdziarski et al. 2008). New campaigns with similar radio resolution, as well as new observations in the ray domain have been obtained since the Dhawan’s et al. original results (Albert et al. 2008b). A key aspect in these highangularresolution campaigns is the observed maintenance in time of the morphology of the radio emission of the system: the changing morphology of the radio emission along the orbit would require a highly unstable jet, which details are not expected to be reproduced orbit after orbit as indicated by current results (Albert et al. 2008b). The absence of accretion signatures in Xrays in Chandra and XMMNewton observations (as reported by Sidoli et al. 2006, Chernyakova et al. 2006, and Paredes et al. 2007) is another relevant aspect of the discussion about the compact object companion. Finally, it is interesting to note that neutrino detection or nondetection with ICECUBE will shed light on the nature of the ray emission irrespective of the system composition (e.g., Aharonian et al. 2006b, Torres and Halzen 2007).
For LS 5039, a periodicity in the ray flux, consistent with the orbital timescale as determined by Casares et al. (2005b), was found with amazing precision (Aharonian et al. 2006). Short timescale variability displayed on top of this periodic behavior, both in flux and spectrum, was also reported. It was found that the parameters of powerlaw fits to the ray data obtained in 0.1 phase binning already displayed significant variability. Current H.E.S.S. observations of LS 5039 ( hours distributed over many orbital cycles, Aharonian et al. 2006) constitute one of the most detailed datasets of highenergy astrophysics. Similarly to LS I +61 303, the discovery of a jetlike radio structure in LS 5039 and the fact of it being the only radio/Xray source colocalized with a mildly variable (Torres et al. 2001a,b) EGRET detection, prompted a microquasar interpretation (advanced already by Paredes et al. 2000). However, the current mentioned findings at radio and VHE rays in the cases of LS I +61 303 (Dhawan et al. 2006, Albert et al. 2006, 2008b) or PSR B125963 (Aharonian et al. 2005), gave the perspective that all three systems are different realizations of the same scenario: a pulsarmassive star binary. Dubus (2006a,b) has studied these similarities. He provided simulations of the extended radio emission of LS 5039 showing that the features found in high resolution radio observations could also be interpreted as the result of a pulsar wind. Recently, Ribó et al. (2008) provided VLBA radio observations of LS 5039 with morphological and astrometric information at milliarcsecond scales. They showed that a microquasar scenario cannot easily explain the observed changes in morphology. All these results, together with the assessment of the low Xray state (Martocchia et al. 2005) made the pulsar hypothesis tenable, and the possibility of explaining the H.E.S.S. phenomenology in such a case, an interesting working hypothesis.
High energy emission from pulsar binaries has been subject of study for a long time (just to quote a nonexhaustive list of references note the works of Maraschi and Treves 1981; Protheroe and Stanev 1987, Arons and Tavani 1993, 1994; Moskalenko et al. 1993; Bednarek 1997, Kirk et al. 1999, Ball and Kirk 2000, Romero et al. 2001, Anchordoqui et al. 2003, and others already cited above). LS 5039 has been recently subject of intense theoretical studies (e.g., Bednarek 2006, 2007; which we comment on in more detail below, BoschRamon et al. 2005; Böttcher 2007; Böttcher and Dermer 2005; Dermer and Böttcher 2006; Dubus 2006a,b; Paredes et al. 2006; Khangulyan et al. 2007; Dubus et al. 2007).
In the penultimate paper mentioned in the list above, Khangulyan et al. (2007), and contrary to the assumption here, authors assumed a jet structure perpendicular to the orbital plane of the system. The energy spectrum and lightcurves were computed, accounting for the acceleration efficiency, the location of the accelerator along the jet, the speed of the emitting flow, the inclination angle of the system, as well as specific features related to anisotropic inverse Compton (IC) scattering and pair production. Different magnetic fields, affecting Synchrotron emission, and the losses they produced, were also tested given a large model parameter space. Authors found a good agreement between H.E.S.S. data for some of their models.
In the last of these papers, Dubus et al. (2007) computed the phase dependent lightcurve and spectra expected from inverse Compton interactions from electrons injected close to the compact object, assumed as a likely rotationpowered pulsar. Since the angle at which an observer sees the binary and propagating electrons changes with the orbit (see below), a phase dependence of the spectrum is expected, and anisotropic inverse Compton is needed to compute it. In general, they found that the lightcurve is a good fit to the observations, except at the phases of maximum attenuation where pair cascade emission plays a role. Dubus et al. (2007) do not consider cascading in their models, as we do here. Without cascading, zero flux is expected at a broad phase around periastron, which is not found. This lack of cascading in their model also affects the spectra, which are not reproduced well, particularly at the superior conjunction broad phases of the orbit. They mentioned that both, cascading and/or a change in the slope of the powerlaw injection for the interacting electron distribution could be needed to explain the spectrum in these phases, what we explore in detail in this work.
In order to compute inverse Compton emission from LS 5039, we use, as in previous works, leptons interacting with the star photon field. Geometry is described there with different levels of detail, what influence the results. In general, cascading processes were not taken into account, and the goodness of fitting the H.E.S.S. data is arguable in most cases, both for the lightcurve and spectrum.
In none of the papers mentioned above, the theoretical predictions for the short timescale spectral variability found by H.E.S.S. in 0.1 phase binning was shown and compared with data. We discuss these results from our model below.
In recent Letters (SierpowskaBartosik and Torres 2007, 2008) under the assumption that LS 5039 is composed by a pulsar rotating around an O6.5V star in the day orbit, we presented the results of a leptonic (for a generic hadronic model see Romero et al. 2003) theoretical modeling for the highenergy phenomenology observed by H.E.S.S. These works studied the lightcurve, the spectral orbital variability in both broad orbital phases and in shorter (0.1 phase binning) timescales and have found a complete agreement between H.E.S.S. observations and our predictions. We have also analyzed how this model could be tested by Gammaray Large Area Space Telescope (GLAST), and how much time would be needed for this satellite in order to rule the model out in case theory significantly departs from reality. But many details of implementation which are not only useful for the case of LS 5039 but for all others close massive ray binaries, as well as many interesting results concerning the binary geometry, wind termination, opacities to different processes along the orbit of the system, and further testing at the highest energy ray domain were left without discussion in our previous works. Here, we provide these details, together with benchmark cases that are useful to understand the formation of the very highenergy lightcurve and phase dependent spectra.
The rest of this paper is organized as follows. Next Section introduces the model concept and its main properties. It provides a discussion of geometry, wind termination, and opacities along the orbit of the system (we focus on LS 5039). An accompanying Appendix provides mathematical derivations of the formulae used and useful intermediate results that are key for the model, but too cumbersome to include them as part of the main text. It also deals with numerical implementation, and describes in detail the Monte Carlo simulation of the cascading processes. The results follow: Section 3 deals with a monoenergetic interacting particle population, and Section 4, with powerlaw primary distributions. Comparison with H.E.S.S. results is made in these Sections and details about additional tests are given. Final concluding remarks are provided at the end.
2 Description of the model and its implementation
Under the assumption that the pulsar in the binary is energetic enough to prevent matter from the massive companion from accreting, a termination shock is created in the interaction region of the pulsar and donor star winds. This is represented in Fig. 1. We focus on the specific case of the binary LS 5039, which we use as a testbed all along this paper. The volume of the system is separated by the termination shock, which structure depends on features of the colliding winds: it may be influenced by the anisotropy of the winds themselves, the motion of the pulsar along the orbit, turbulences in the shock flow, etc. For simplicity it is assumed here that the winds are radial and spherically symmetric, and that the termination of the pulsar wind is an axial symmetric structure with negligible thickness. In this general picture there are three regions of different properties in the binary: the pulsar wind zone (PWZ), the shock (SR), and the massive star wind zone (MSWZ).
The energy content in the interaction population of particles is assumed as a fraction of pulsar spindown power . In case of young energetic pulsars, this power is typically . Propagating pairs upscatter thermal photons from the massive star due to inverse Compton process. For close binaries, the radiation field of hot massive stars (type O, Be or WR, having typical surface temperatures in the range and linear dimension ) dominates along the whole orbit over other possible fields (e.g., the magnetic field or the thermal field of the neutron star). This thermal radiation field is anisotropic, particularly for injected close to the pulsar (the radiation source is misplaced with respect to the electron injection place).
The highenergy photons produced by pairs can initiate cascades due to subsequent pair production in absorption () process with the same radiation field (as sketched in Fig. 1). We assume that these cascades develop along the primary injection direction, i.e., in a onedimensional way, which is certainly justified based on the relativistic velocity of the interacting electrons. This process is followed up to the termination shock unless leptons lose their energy before reaching it. Those leptons which propagate to the shock region are trapped there by its magnetic field. Radiation from them is isotropised. The photons produced in cascades which reach the shock can get through it and finally escape from the binary or be absorbed in the radiation field close to the massive star, some may even reach the stellar surface.
In the shock region leptons move along it with velocity . They could be reaccelerated and produce radiation via synchrotron (local magnetic field from the pulsar side) or inverse Compton scattering (ICS, thermal radiation field from the massive star) processes. However, as they are isotropised in the local magnetic field, photons are produced in different directions and their directionality towards the observer is lost. It was already shown by Sierpowska and Bednarek (2005) that in compact binary systems (as an example, the parameters of Cyg X3 were taken by these authors) the radiation processes in the shock region do not dominate: the energy carried by reaching the shock is a small fraction of total injected power. Furthermore, we will show that for the parameters relevant to the LS 5039 scenario, the PWZ is relatively large with respect to the whole volume of the system for the significant range of the binary orbit.
2.1 Hydrodynamic balance
Assuming that both, the pulsar and the massive star winds are spherically symmetric, and based on the hydrodynamic equilibrium of the flows, the geometry of the termination shock is described by parameter where and are the loss mass rates and velocities of the two winds (Girard and Wilson, 1987). The shock will be symmetric with respect to the line joining two stars, with a shock front at a distance from the one of the stars:
(1) 
The surface of the shock front can be approximated then by a conelike structure with opening angle given by where takes the smaller value between the two magnitudes quoted. This last expression was achieved under the assumption of nonrelativistic winds in the simulations of the termination shock structure in Girard and Wilson (1987), albeit it is also in agreement with the relativistic winds case (e.g., Eichler and Usov, 1993; Bogovalov et al., 2007). If one of the stars is a pulsar of a spin down luminosity and the power of the massive star is , the parameter can be calculated from the formula (e.g., Ball and Kirk 2000). Note that for , the star wind dominates over the pulsar’s and the termination shock wraps around it. Note also that for , the shock is at equal distance, , between the stars. In the case of LS 5039, for the assumed (nominal) spindown luminosity discussed below, the value of is between (periastron) and (apastron).
The massive star in the LS 5039 binary system is of O type, which wind is radiation driven. The velocity of the wind at a certain distance can be described by classical velocity law where is the wind velocity at infinity, the hydrostatic radius of the star, is the velocity close to the stellar surface, and (e.g., Cassinelli 1979, Lamers and Cassinelli 1999) and we assume . As can be seen by plotting the velocity law, the influence of this parameter is minor. Typical wind velocities for O/Be type stars are . For LS 5039 we have and (Casares et al. 2005). The assumed value for the star massloss rate is M yr, while the typical values for O/Be stars are M yr.
2.2 The PWZ and the interacting lepton population
The magnetization parameter , (e.g., Langdon et al. 1988) is defined as the ratio of Poynting flux to relativistic particle energy flux. The magnetic field in is that of the upstream shock propagating with bulk Lorentz factor and being the relativistic particle density. The processes establishing the effective change of along the PWZ are the central issue in the discussion of dissipation mechanisms in relativistic plasma flows, exemplified with the Crab pulsar wind (e.g., Kennel and Coroniti, 1984), where variations seem notable. The Crab wind is originally Poyntingdominated ( close to the neutron star, e.g. Arons 1979); but it is kineticdominated near the termination shock (, e.g., Kennel & Coroniti 1984). The change in is produced as a result of dissipative plasma processes in the PWZ, which is characterized by high bulk Lorentz factor (e.g., Melatos 1998). Dissipation (plasma processes engaged in the conversion of electromagnetic towards particle kinetic energy) in Poynting flux dominated plasma flows can be in the form of stochastic/nonstochastic and adiabatic/nonadiabatic processes, thermal heating/non thermal particle generation, and isotropic adiabatic expansion/directed bulk acceleration of the plasma flow (Jaroschek et al. 2008). The microphysical details are decisive when looking for the type of particle energization and their spectra.
The existence of wisps in the inner structure of the Crab nebula has been discussed by, e.g., Lou (1998). The interesting fact of some of them being close to the pulsar, apparently well inside the PWZ, was interpreted as being produced by slightly inhomogeneous wind streams, demonstrating that reverse fast MHD shocks at various spin latitudes can appear quasistationary in space when their propagation speeds relative to the pulsar wind are comparable to the relativistic outflow. The possibility of an inhomogeneous wind stream is not implausible. Successive radio pulses from a pulsar indeed vary in their shapes, eventhough the average pulse is stable. A slower wind stream will be eventually caught by a faster one to trigger forward and reverse fast MHD shocks inside the PWZ. In this zone, charged particles can be further accelerated by these turbulences and magnetic reconnection. Lou proposed that an isotropised powerlaw like energy distribution of the electrons thus produced help to understand the properties of the changing and brilliant inner nebula.
A monoenergetic assumption for the distribution of leptons in the PWZ can be considered as a first approach to the problem. On one hand, the magnetization parameter may be a function of angle, and although must be very small in the equatorial part of the wind (e.g., Kirk 2006), simulations do not favor an angleindependent low value (Komissarov & Lyubarsky 2004). On the other hand, Contopoulos & Kazanas (2002) already showed that the Lorentz factor of the outflowing plasma could increase linearly with distance from the light cylinder (implying that decreases inversely proportional to the distance). Contopoulos & Kazanas (2002) mentioned that this specific radial dependence of the pulsar winds Lorentz factor is expected to have additional observational consequences: e.g., Bogovalov & Aharonian (2000) computed the Comptonization of soft photons to TeV energies in the Crab through their interaction with the expanding MHD wind, while Tavani & Arons (1997) and Ball & Kirk (2000) computed the corresponding radiation expected by the radiopulsar Be star binary system PSR B125963 through the interaction of the relativistic wind with the photon field of the companion in much the same way we do here. Indeed, the details in these predictions would be modified, as shown by Sierpowska & Bednarek (2004, 2005) should the linear acceleration model be adopted. Hibschman and Arons (2001) discussed the creation of electronpositron cascades in the context of pulsar polar cap acceleration models. They computed the spectrum of pairs that would be produced outflowing the magnetosphere. They found that the pair spectra should be described by a powerlaw.
One possibility for the dissipative conversion is established by magnetic reconnection processes between antiparallel magnetic stripes during outwards propagation in the PWZ (Lyubarsky and Kirk, 2001; Kirk and Skjaeraasen, 2003). Kirk (2004) considered acceleration in relativistic current sheets (large magnetization parameter, with Alfven speed close to c). Recently, Jaroschek et al. (2008, and see references therein for related work) addressed the problem of interacting relativistic current sheets in selfconsistent kinetic plasma simulations, identifying the generation of nonthermal particles and formation of a stable powerlaw shape in the particle energy distributions . Depending on the dimension of the simulation, spectral index from 2 (1D, attributed to a stochastic Fermitype acceleration) to 34 (recognized as a rather universal index of relativistic magnetic reconnection in previous 2D and 3D kinetic simulations, see Jaroschek et al. 2004, Zenitani and Hoshino 2005) were found. Lyubarsky and Liverts (2008) also studied the compression driven magnetic reconnection in the relativistic pair plasma, using 2.5D (i.e., 2D spatial, 3D velocity) simulations, finding that the spectrum of particles was nonthermal, and a power law was produced. It seems a powerlaw distribution for the leptons inside the PWZ is then a plausible assumption.
All in all, to find an apriori dissipation solution for pulsars, and in particular, for the assumed pulsar in the LS 5039 which is the one we focus, is beyond the scope of this work (and actually, for the latter particular case, such solution is beyond what is by definition possible for a pulsar that we do not know exists). In general, we note that an additional difficulty resides in the fact that the PWZ of pulsars in binaries may be subject to conditions others than those found in isolated pulsars. It is not implausible that close systems may trigger different phenomenology within the PWZ, ultimately affecting particle acceleration there. Nevertheless, it is relevant for this paper to assume a particular particle’s energy spectra with which we compute high energy processes in the PWZ, e.g. the upComptonization of the stellar field. We will assume two cases, a monoenergetic spectra as a benchmark (e.g., see Bogovalov & Aharonian 2000) and a generic powerlaw spectra (that could itself be subject to orbital variability). Both act as a phenomenological assumption in this paper, which goodness is to be assessed a posteriori, by comparison with data.
As discussed, initial injection could come directly from the pulsar (the interacting particle population can of course be later affected by the equilibrium between this injected distribution and the losses to which it is subject, just as in the case of shockprovided electron primaries). The more compact the binary is, the more these two settings (shock and pulsar injection, equilibrated by losses) are similar to each other. Given the directionality of the high Lorentz factor inverse Compton process, photons directed towards the observer are generally those coming from electrons moving in the same direction. Opacities to processes such as inverse Compton and absorption are high in close ray binaries, cascades can develop, and highenergy processes can already happen, as we explicitly show below, in the preshock region.
In the model where pairs are injected as monoenergetic particles with the energy corresponding to the bulk Lorentz factor of the pulsar wind, they are frozen in the Bfield. Under this assumption we neglect here the synchrotron losses, since there are none.
When the injected pairs distribution is given by a power law spectrum the situation is more complex, since not all of them may be frozen in the PWZ field. The synchrotron cooling time is given by equation:
(2) 
where is the local magnetic field and is given in Gauss, and is electron energy given in TeV. For IC scattering, the timescale is instead given by
(3) 
which gives good approximation for electron energy loss time for TeV (e.g., Khangulyan et al. 2008). Timescales approach when
(4) 
For a star with effective temperature K, the thermal field density at certain point at distance from the massive star center is given by:
(5) 
Thus, the local magnetic field in the pulsar wind region is given by:
(6) 
Applying this condition at periastron (), the local magnetic field at the assumed injection place results in G. At apastron (), it results in an stronger condition G. That is, the magnetic field should be less than these values in order for IC to dominate over synchrotron losses at that particular position (the light cylinder) in the PWZ. Figure 2 shows this in detail.
The magnetic field at the light cylinder distance is given by the dipole formula . In the PWZ, the magnetic field is decreasing with distance as . For a millisecond pulsar we get cm and G (assuming G and ms) up to cm and G (assuming G and ms), where km. The results for these different parameters are similar as there were obtained with fixed pulsar power in the model erg s and they are related by the standard formula . We also assume here that the magnetization parameter is , but have explored other values of this and other parameters as well, with similar results (see Figure 2).
Given our results (see Fig. 2) where we show the local magnetic field for which IC dominates and the magnetic field in the PWZ as a function of the distance from the light cylinder we can conclude that the injection for the model have to (generically) occur at some distance from the light cylinder, or/and, if closer to it, the synchrotron losses can be important. However, as the separation of the binary is cm, several orders of magnitude larger than the light cylinder (e.g., cm), the change of the injection place within e.g. already gives the initial injection distance at . Thus, no significant effect in the PWZ photon spectra and lightcurve is produced.
2.3 Normalization of the relativistic particle power
The fraction of the pulsar spindown power ending in the interacting pairs can then be written as:
(7) 
Assuming that the distance to the source is kpc, the normalization factor for electrons traveling towards Earth is The specific normalization factors in the expression of the injection rate will be given together with the results for two models in the corresponding sections below. In the models presented here, only a small fraction (1%) of the pulsar’s ends up in relativistic leptons. This is consistent with ions carrying much of the wind energetics. In the case of monoenergetic lepton distribution, where the energy of the primaries is fixed at TeV, we have . In the case of a powerlaw in energy, that may be constant or vary along the orbit, we have .
2.4 On parameter interdependencies
The expression represents the realization of a specific scenario fixed by different values of the pulsar and massive star parameters. For instance, assuming and as given in the Table 4 below, we get for the periastron phase. But we would achieve the same value of with a smaller , say 10 , if within uncertainties we adopt a smaller value of . Clearly, if only one of the quantities, or changes, the shock position moves, what results in a correspondingly smaller or larger PWZ in which cascading processes are set. A key parameter is then given by , instead of , since this is what defines the real size of the PWZ. The mild dependence of on makes possible to accommodate a large variation in the latter maintaining similar results. Fig. 3 shows this dependence for a relevant range of .
In addition, it is obvious that and are linearly (inversely) related. We note that a fixed value of does not imply a specific value for as it is combined with parameters of the massive star wind, subject to uncertainty in their measurements for the specific case of binary treated (e.g., for LS 5039), if at all known. In the end, the same results can be obtained for different sets of parameters. There is a degeneracy between the shock position, which determines the extent of the PWZ, and the injected power in relativistic electrons. For a smaller PWZ (e.g., if we assume a smaller for the same product of ) a larger amount of injected power compensates the reduced interacting region. We found that to get the similar results when the parameter is smaller, the parameter have to be increased roughly by the same factor.
2.5 Wind termination
Very interesting in the context of LS 5039 model properties is the dependence of the distance from the pulsar to the termination shock in the direction to the observer as a function of the orbital phase. This is shown in Fig. 4. We find that for both inclinations, the wind is unterminated for a specific range of phases along the orbit, i.e., the electron propagating in the direction of the observer would find no shock. The PWZ would always be limited in the observer’s direction only if the inclination of the system is close to zero, i.e., the smaller the binary inclination the narrower the region of the wind nontermination viewed by the observer. In the discussed scenarios for LS 5039, the termination of the pulsar wind is limited to the phases for inclination and to for , what coincides with the phases where the emission maximum is detected in the very highenergy photon range (see below the observed lightcurve obtained by H.E.S.S.). The unterminated wind is viewed by the observer at the phase range between apastron and INFC, while a strongly limited wind appears from periastron and SUPC. Note that the important differences in the range of phases in which the unterminated wind appears for distinct inclinations (e.g., the observer begins to see the unterminated pulsar wind at for compared to for ) produce a distinguishable feature of any model with a fixed orbital inclination angle.
From the geometrical properties of the shock surface, there can also be specific phases for which the termination shock is directed edgeon to the observer. These are phases close to the nonterminated wind viewing conditions: slightly before and after those specific phases for which the wind becomes nonterminated, e.g., for the case of LS 5039 and , and ; whereas for , it appears at and . Thus, we find that the condition for lepton propagation change significantly in a relatively short phase period, when the termination shock is getting further and closer to the pulsar (phases periods and ). In Table 1 we show the differences in these geometrical parameters for a few characteristic orbital phases. Note that even when they are essential to understand the formation of the lightcurve and phase dependent high energy photon spectra, all the above features are based on the simple approximation of the geometry of the termination shock. In a more realistic scenario, the transition between the terminated and unterminated wind, and its connection with orbital phases and inclination are expected to be more complicated yet.
Periastron  0  0.0  2.25  0.8  0.0  2.25  0.6  
SUPC  44  0.06  2.45  0.8  0.06  2.45  0.6  
134  0.27  4.06  3.2  0.27  4.06  3.2  
171  0.45  4.7  …  0.36  4.5  …  
Apastron  180  0.5  4.72  …  0.5  4.72  …  
INFC  224  0.72  4.1  …  0.72  4.1  …  
285  0.89  2.8  …  0.92  2.6  …  
314  0.94  2.47  2.4  0.94  2.47  2.4 
2.6 Opacities along the LS 5039 orbit
The conditions for leptonic processes for this specific binary can be discussed based on optical depths to ICS and absorption. The target photons for IC scattering of injected electrons and for pair production for secondary photons are low energy photons of blackbody spectrum with temperature , which is the surface temperature of the massive star. This is an anisotropic field as the source of thermal photons differ from the place of injection of relativistic electrons, which for simplicity is assumed to be at the pulsar location for any given orbital phase. For fixed geometry parameters, the optical depths change with the separation of the binary (in general with the distance to the massive star), the angle of injection (the direction of propagation with respect to the massive star), and the energy of injected particle (the electrons or photons for absorption).
To have a first handle on opacities, we have calculated the optical depths adopting the binary separation at periastron and apastron. In case of the optical depths for photons we had to do an additional simplifying assumption to allow for a direct comparison with the optical depth for electrons. In the discussed model of this paper, photons are secondary particles and they do not have the same place of injection as electrons, but because the path of electrons up to their interaction is much smaller than the system separation, we assume now the same place of injection for photons and electrons: for a first generation of photons, this is sufficient for comparison.
Results are shown in Fig. 5 and the needed formulae for its computation are given in the Appendix.
The optical depths for ICS can be as high as few for significant fractions of the orbit, and are still above unity for energies close to TeV. With respect to the injection angle, the opacities are above for all propagation angles, except the outward directions at the apastron separation. (Recall that outward –inward– directions
refer to angles close to 0 and 180, respectively, as shown in Fig. 1).
In the case of
pair production, the interactions are limited to specific energies of the photons and are strongly angledependent. The most favorable case for pair creation is represented by those photons with energy in the range to few TeV, propagating at .
The highest optical depths to pair production are found for directions tangent to the massive star surface: for propagation directly towards the massive star pair production processes are limited by the star surface itself (for periastron the tangent angle is at ; for apastron, it is given at an angle of ). This is better shown in Fig. 6.
The opacities for pair production are characterized by maxima which change with respect to the angle of photon injection. For inward directions (in general for angles larger then ), the peak energy is higher, , than for outward directions, where we find . For energies up to the peak for pair production, the IC process dominate for the same angle of particle injection, but for higher energies the opacities for absorption become slightly higher than for ICS.
Cascading processes are effective when the interaction path for electrons and photons are short enough so that a cascade can be initiated already inside the pulsar wind zone. As we can see in Fig. 5, the opacities at highenergy range (KleinNishina) are comparable to the maximal opacities for absorption. This means similar probability of interaction for both photons and electrons. If the electrons scatter in the KleinNishina regime, the photon produced, with comparable energy to that carried by the initiating electron, is likely to be absorbed, and a cascade is produced.
The dependence of the optical depths upon the orbital phase for specific parameters of the LS 5039 is shown in Fig. 7. These opacities are calculated up to the termination shock in the observer direction. The presence of the shock (at a distance ) limits the optical depths. Since this parameter is highly variable along the orbit (see Fig. 4), its influence on the optical depth values is not minor. As the angle to the observer, which defines the primary injection angle, vary in the range (for INFC and SUPC, respectively), we found that there is a range of orbital phases for which the PWZ is nonterminated. In that case the cascading process –which develops linearly– is followed up to electron’s complete cooling (defined by the energy ). The optical depths in the nonterminated pulsar wind are thus comparable to those presented in Fig. 5, taking into account the differences between the separation and angle of injection.
Apart from the condition for electron cooling, in the case of the terminated wind the cascade is followed up to the moment when the electron reaches the shock region, whatever happens first. As can be seen from comparison of Fig. 5 and 7, the opacities are significantly smaller when the propagation is limited. For instance, as a comparison we can choose the electron injection at periastron, the angle ( for ), and fix the energy to TeV. Then, the opacities along the orbit up to the terminated shock are while for the unterminated processes the opacity are . A direct comparison for specific energies of injected particles and the two inclinations considered is shown in Fig. 8. Moreover, the opacities for both processes decreases along the propagation path, see Fig. 9, as the angle to the observer also decreases. This indicate that even if the photons propagate very close to the massive star, the probability for absorption can be much smaller than that at its injection place, and finally become less then at the distance from the pulsar , for an energy of TeV and TeV respectively.
3 Monoenergetic electrons in the PWZ
Following the normalization formula given by Eq. (7), for the primary lepton energy TeV and the nominal value of pulsar spindown power , and assuming a distance to the source kpc we get The normalization factor for electrons traveling towards Earth is . As a first normalization of the simulation results below, we have fixed .
3.1 Lightcurve and broadphase spectra
Meaning  Symbol  Adopted value 

Spindown power of assumed pulsar  erg s  
Fraction of in relativistic electrons  
The energy of injected pairs  10 TeV 
Based on the simulations for the specific phases along the orbit, the photon spectra and lightcurves in different energy ranges ( GeV, GeV, and above TeV) were calculated for two assumed inclinations of the orbit of the system, and . The highest energy range corresponds to the data presented by H.E.S.S. (Aharonian et. al. 2006). To allow for a comparison, the photon spectra were calculated also for a primary electron energy of TeV for specific orbital phases (periastron, apastron, SUPC and INFC).
The photon spectra (SED) are presented in Fig. 10. All the spectra for primary electron energy TeV are hard, presenting a photon index , and they have significantly higher flux at the highest energy range, close to the initial primaries’ TeV, as the ICS of primary pairs mainly occur in KleinNishina range. Below this peak the spectra are well represent by the powerlaw with photon index close to , albeit one can also notice the change in the photon index from SUPC to INFC (harder spectrum) and the notably different energy cutoffs. These features are discussed below.
While for phases around SUPC the PWZ is limited in the direction to the observer, for the opposite phases the PWZ is unterminated. On the other hand, for INFC the angle to the observer is at its minimum (and also the separation of the system is larger) what causes the decreasing of optical depth to both processes considered. In contrary to the INFC phases, those at SUPC present strong absorption features at an energy range from to few TeV, which cause the spectra to be cut at lower energies. This is due to the dependence of the optical depth to absorption on the energy of photon and the influence of geometry in defining the photon path towards the observer. The larger the angle to the observer, the stronger the effects of absorption are, what can be seen from comparison of the results for two inclination angles (SUPC and periastron phases). Notice, that for the angle to the observer attains the largest value of the discussed examples, . The absorption for SUPC phases takes place already in the PWZ, what can initiate the cascading processes. Photons which pass through the shock region can be also absorbed in the MSWZ. The opacities for absorption also strongly depend on photon energy (see Fig. 8). High energy photons, with energy from few GeV up to few TeV, are the most likely to be absorbed, as the opacities decrease with energy for all shown angles of the photon injection. Indeed, the number of the photons absorbed in MSWZ is about % of the number of photons reaching the shock region and the effect is significant in the final spectra. On the other hand, the cascading in PWZ produce lower energy leptons and photons (in the cooling process of secondary pairs) and this causes a higher photon flux at energies up to few 10 GeV. The processes in PWZ are limited by the presence of the shock, what influence the photon production rate. For INFC phases the absorption of rays is minor, as the photons propagate outwards of the massive star; once produced, most photons can escape from the system (also because of the threshold to pair production). Higher energy photons are produced mostly in the first IC interaction of the primary pairs, while further electron cooling, not limited by the termination shock, supply the spectrum in lower energy photons. This flux is not as high as in the case of SUPC phases, where most of the cascading takes place.
Similar dependencies, both for INFC and SUPC phases, are present in the spectra obtained from the simulations for primary energy of pairs TeV. As the optical depths to IC scattering are higher in that case, the processes are more efficient and the number of produced photons are higher. To give an example, most photons of energy 2030 GeV at SUPCperiastron phase are absorbed in MSWZ.
Comparing the photon spectra produced at SUPC and INFC, an anticorrelation between GeV and TeV photon fluxes is evident, at least comparing the spectra at energies below and above GeV. All these effects play a role in the formation of the ray lightcurve, shown in Fig. 11. The lightcurve in the highest energy range ( TeV) has a broad minimum around SUPC, . From the comparison of the position of the termination shock along the orbit, the opacities to ICS and absorption we see that this minimum is formed despite being at phases with the highest opacities to both processes considered and so having an effective cascading in the PWZ region. This minimum is mainly due to the absorption of the photons which get through the shock and propagate into the massive star wind region. The broad maximum in the highenergy lightcurve is formed at opposite phases around INFC, . The maximum for higher inclination is also characterized by a local minimum close to the INFC phase, what is the result of the IC opacity dependence on the propagation angle. The range of INFC phases corresponds to the unterminated pulsar wind in the observer direction. The opacities for both inclinations have local peaks in this range, what was discussed in the previous paragraph. The local minimum in photon flux within the INFC range reflects a similar behavior of the opacities. So, as far as the propagation of the particles in the PWZ is unterminated, the highenergy lightcurve formation is in agreement with the dependence of the optical depths. Additionally, we can also see that the first local peak in this lightcurve is formed earlier in phase for inclination (), what is also in agreement with the dependence shown by the distance to the termination shock. Note that the second peak in the broad TeV lightcurve maximum, at , is higher than the first one, at what is the result of a change in the separation of the system together with the angle to the observer.
3.2 Peaks and dips in the lightcurve
The geometrical conditions for leptonic processes change significantly along the orbit, and they are more efficient for phases around SUPC. However, these processes, as we can already see from the opacities dependence on the orbital phases, is limited by the presence of the termination shock. This causes that the best conditions for very highenergy photon production, and finally escaping from the system, occur for an specific combination of the angle to the observer, the separation of the system, and the distance to the shock from the pulsar side. From Figs. 7 and 4 we can see that this happens at the phases and what reflects in the TeV lightcurve. At phases close to INFC photons are produced mainly in the primary pairs cooling when the propagating electron undergo frequent scatterings. Together with the fact that there is no efficient absorption of photons once they are produced, they finally escape from the system, what yields to the broad maximum in the TeV lightcurve. On the other hand, at phases around SUPC, highenergy photons are absorbed when propagating through the system in the MSWZ, what causes a dip in the lightcurve. For these phases many more lower energy photons are produced in cascades in the PWZ, what yields to a maximum in the GeV lightcurve for SUPC, anticorrelated with the behavior at TeV energies.
4 A powerlaw electron distribution
For further exploration of the ray production model we set a new assumption: the energy distribution of the interacting pairs is given by a powerlaw spectrum. We will additionally assume that the powerlaw may be constant or vary along the orbit. Motivated by the different observational behavior found, we have assumed that two different spectral indices correspond to the two broad orbital intervals proposed by HESS (Aharonian et. al., 2006). For direct comparison we have specified the interval around the inferior conjunction (with the apastron phase): , and around superior conjunction (including the periastron phase): and as being bathed by different electron distributions. The results for powerlaws distribution were already summarily presented in our earlier work SierpowskaBartosik and Torres (2007a,b) and we refer to these works for further details. The agreement in both spectra and lightcurve is notable (particularly for the case of a variable lepton spectrum along the orbit), as can be seen in Figs. 14 and 15 discussed below.
Meaning  Symbol  Adopted value 
Spindown power of assumed pulsar  erg s  
VHE cutoff of the injection spectra  TeV  
Constant lepton spectrum along the orbit  
Fraction of in leptons  
Slope of the powerlaw  
Variable lepton spectrum along the orbit  
Fraction of in leptons at INFC interval  
Slope of the powerlaw at INFC interval  
Fraction of in leptons at SUPC interval  
Slope of the powerlaw at SUPC interval 
As we could have already noticed from Fig. 5, the optical depth for high energy electrons is below unity. In that case, part of the initial electrons will be interacting in the PWZ less efficiently and finally will reach the shock region (we remind that in the direction of the observer, there is not always a shock in the electron’s propagation). In Fig. 12 we show the spectra of initial electron distribution and corresponding electron spectra after propagation in the PWZ. The electrons which reach the termination shock are isotropised there. The termination shock shape is specific for each phase and it is limited in space. The electrons at the shock, locally reaccelerated, become the initial spectra for the next generation of photons (not only radiative but adiabatic losses have to be taken into account).
The MSWZ can play a role (although we believe it will not be too important in such close binaries like LS 5039, because of the high opacities already encountered in the PWZ along most of the orbit and the loss of directionality –consequently, of random photon emission– of the electron population). An estimation of the contribution of the cascades in MSW is not trivial, especially from the normalization point of view, as the magnetic filed in this region causes isotropisation of produced photons. In addition, the magnetic field of the massive star is ordered, what could cause additional effects e.g. focusing the propagating electrons in some regions close to the massive star (where the magnetic field is dipolar) as was shown in Sierpowska & Bednarek (2005). The impact of the MSWZ could make the need for a change in the injection index of relativistic particles less severe. A 3D cascading code is needed for such estimates, and we expect to report on that in the future.
The model assumes that the initial electron are injected in the vicinity of the pulsar light cylinder. This is actually an assumption which we have investigated further. Indeed. we have investigate also the changes in the produced photon spectra if injection take place at a further distance inside the PWZ. Results are given in Fig. 13. For phases where the wind is terminated in the direction to the observer (e.g., SUPC) the new injected radius was fixed to from the massive star, where is separation at given phase (i.e., at the middle of the PWZ). For phases such as INFC, where the shock is unterminated, the injection place was shifted by the radius of the star , that also corresponds to the halfdistance of the separation between the pulsar and shock for this phase, if movement is in the direction to the massive star. We can see that the produced photon spectra do not differ significantly in this two scenarios, making our results stable.
4.1 Lightcurve and broadphase spectra
To construct the lightcurve, spectra for over 20 orbital phases were calculated to cover the whole orbit of LS 5039. The specific averaged spectra were obtained based on the same orbital intervals as presented in H.E.S.S. data. The averaged spectra were obtained summing up individual contributions from orbital spectra in given phases, each with a weight () corresponding to the fraction of orbital time that the system spends in the corresponding phase bin. Then, comparing the observational and theoretical averaged spectra, the parameter was estimated. Both for constant and variable injection model, we get that the fraction of the spindown power in the primary leptons has to be at the level of %. With this parameter in hand each of the single spectra can be equally normalized.
It is worth noticing how well these lightcurves compare with those in the work by Bednarek (2007), at least for some of the specific phases considered by him. Bednarek also included cascading in his simulations, and the physical input of his model (although in the case of a microquasar scenario) is similar to ours. As a result, the anticorrelation phenomena (from GeV to TeV energies) is also a result of his work. The spectrum along the orbit with respect to H.E.S.S. datapoints and the possible shorttimescale variability (see below) was not provided by Bednarek, so that a comparison with these results is not possible. A detailed differentiation between these two models is a key input for distinguishing (microquasar or ray binary) scenarios, even when some assumptions are intrinsic to each of the models, isolating the contribution of an equal physical input can help decide on what object constitute the system LS 5039. As an example: Bednarek did find in his model that a fixed inclination angle ( = 60) was needed in order to reproduce the shape of the H.E.S.S. lightcurve results, whereas in our case, as we see in Figure 15 the influence of inclination is minor.
For completeness, we mention that as noted above, H.E.S.S. has also provided the evolution of the normalization and slope of a powerlaw fit to the 0.2–5 TeV data in 0.1 phasebinning along the orbit (Aharonian et al. 2006). The use of a powerlaw fit was limited by low statistics in such shorter suborbital intervals, i.e., higherorder functional fittings such as a powerlaw with exponential cutoff were reported to provide a no better fit and were not justified. To directly compare with these results, we have applied the same approach to treat the model predictions, i.e., we fit a powerlaw in the same energy range and phase binning. We show here this comparison in Fig. 16 (taken from SierpowskaBartosik and Torres 2008) in the case of a variable lepton distribution along the orbit. We find a rather good agreement between model predictions and data. Results for constant lepton distribution along the orbit can be seen in SierpowskaBartosik and Torres (2007a).
In Fig. 14 the SED in H.E.S.S. energy range are shown for both, constant and variable lepton distributions along the orbit, and two inclination of the system and . It was shown in the previous Section that in case of the monoenergetic injection, there are photon flux differences in INFC phase due to the change in the angle to the observer. In the powerlaw distribution model this effect is less significant due to the contributions of different energies, and the dependence on energy of the processes involved. Based for instance on the injection of constant spectrum with slope we can see that there is no significant differences involving the inclination.
Exploring the more evolved model of the variable injection we see that the steeper the primary spectrum is, the higher the photon flux at lower energy range (SED for lower energies were shown by SierpowskaBartosik and Torres 2008). For constant lepton distribution along the orbit, the differences in photon flux below 100 GeV are not so significant, but the difference gets important in the case of a variable lepton distribution (the difference between both models in the lower energy range is about one order of magnitude). Overall, the variable lepton distribution provides a better agreement with all data (note that it has only two extra free parameters when compared with the constant distribution case, see Table 3, but matches more than 10 data points that were missed in the previous case). In particular, it is worth noticing that there is good agreement with the H.E.S.S. spectra at an energy GeV where both SUPC and INFC spectra coincide and above which the photon flux for INFC is higher then the flux for SUPC dominating in the lower energy range; i.e., the energy where the anticorrelation begins (see SierpowskaBartosik and Torres 2008). It is also interesting to remind that absorption alone would produce strict modulation (zero flux) in the energy range 0.2 to 2 TeV whereas the observations show that the flux at TeV is stable. Additional processes, i.e. cascading, must be considered to explain the spectral modulation. Our model have these processes consistently included and the lightcurve details arise then as the interplay of the absorption of rays with the cascading process, in the framework of a varying geometry along the orbit of the system.
4.2 Testing with future data
Apart from possible testing with GLAST (at the level of lightcurve, spectra, hardness ratios, and differentiation between constant and variable electron distribution, see SierpowskaBartosik and Torres 2008), we can provide further possible tests at high ray energies. It was already said that powerlaws do not always present the best fit to the specific spectra along the orbit. Even when fitting such powerlaws to the theoretical predictions provides agreement with data (see Fig. 16), observations with larger statistics (with H.E.S.S., H.E.S.S. II, or CTA) could directly test the model in specific phases. A model failure in specific phases would allow further illumination about the physics of the system. Figure 17 shows the evolution of individual spectrum in the best fitting case of variable lepton distribution along the orbit, at individual phase bins, from 0.1 to 0.9, for testing with future quality of data.
5 Concluding remarks
We presented the details of a theoretical model for the highenergy emission from close ray binaries, and applied it to the particular case of LS 5039. The model assumes a pulsar scenario, where either the pulsar or a closetothepulsar shock injects leptons that after being reprocessed by losses to constitute a steady population, are assumed to interact with the target photon field provided by the companion star within the PWZ. The model accounts for the highly variable system geometry with respect to the observer, and radiative processes; essentially, anisotropic KleinNishina ICS and absorption, put together with a Monte Carlo computation of cascading. The formation of lightcurve and spectra in this model was discussed in detail for the case where the interacting leptons are assumed monoenergetic and described by powerlaws.
Comparing the interacting models of monoenergetic leptons and powerlaw distributions we can see similar dependencies for the spectral changes along the orbit and the GeV to TeV lightcurves. Thus, the case of a monoenergetic population is useful to understand some of the aspects regarding the formation of the observational features, although it does not match observational data. For the monoenergetic case, we have shown the spectra produced at characteristic phases (INFC, SUPC, periastron, apastron) for primary energies of 1 TeV and 10 TeV. In the powerlaw distribution model, the specific features of the photon spectra and the lightcurves produced for an specific primary electron energy overlap. We can still notice the absorption features in the spectra produced at SUPC. We have discussed effects solely based on the optical depths and the general geometrical dependence along the orbit, as well as on the presence of the shock which terminates the pulsar wind in the direction to the observer at SUPC phases.
A powerlaw lepton distribution interacting in the PWZ describes very well the phenomenology found in the LS 5039 system at all timescales, both flux and spectrumwise, even at the shortest timescales measured. This latter result is unexpected: we find that there is nothing a priori in the model that allows one to predict that when broad phase spectra data (INFC and SUPC) are reproduced so will be the data at the individual and much shorter phasebinning, less with such a good agreement. This result point perhaps to some reliability of the model, at least in its essential ingredients: geometry, cascading, interacting electron population.
However, this model certainly has room for improvement. We emphasize here that we do not have an a priori model for the interacting lepton population itself, although we have discussed the research on dissipation processes in the PWZ which may give raise to such distributions if it results from pulsar injection. In any case, the assumption of powerlaws is an approximation to a more complex scenario where the real interacting lepton population is the result of a full escapeloss equation.
In addition, we are not considering yet the multiwavelength emission at lower energies, since we left out of our description the synchrotron emission of electrons accelerated at the shock and the morphology of the shock along the orbit, what we expect to discuss elsewhere.
Other than system scalings that are fixed by multiwavelength observations, the model is based on just a handful of free parameters, and it is subject to tests at high and very highenergy ray observations with both GLAST (described in more detail in SierpowskaBartosik & Torres 2008) and future samples of data at higher energies, where more statistics at finer phase bins can determine better the spectral evolution along the orbit of this interesting system.
We acknowledge extended use of IEECCSIC parallel computers cluster. We acknowledge W. Bednarek for discussions, and the Referee for useful comments. This work was supported by grants AYA 200600530 and CSICPIE 200750I029.
Appendix: Numerical implementation and formulae
This Appendix introduces further essential details concerning geometry and formulae for the implemented process of Inverse Compton scattering and pair production in the anisotropic radiation field of the massive star.
MonteCarlo implementation for the cascading process
A MonteCarlo procedure is applied to calculate the place of electron interaction, , and the energy of the resulting photon in the ICS, ; as well as the place of photon interaction, , and the energy of the produced electron/positron, (see, e.g., Bednarek 1997). The following summarizes the procedure.
When the primary electron is injected, the computational procedure for IC scattering is invoked first in order to get the initial place of interaction (the radial distance from the injection position) and the energy of the upscattered photon if the interaction takes place. For this same electron, the procedure is repeated afterwards up to the moment of electron cooling or when reaching the shock region. The photons produced along the electron path switch on a parallel computational procedure for pair production. This in turn gives the place of photon absorption and the energy of the pair if the process occurs. Because the photons can get trough the termination shock if not absorbed in the PWZ, the numerical procedure is not limited to the termination shock radius. If a pair is created within the PWZ, then the procedure for IC scattering is initiated for it, from which a next generation of photons can be produced. When photons cross the shock, information about them is separately saved in order to get account of the level of flux absorption in the MSWZ. The emission spectra are produced from photons which finally escape from the binary, so the photons absorbed in the MSWZ are not included in it.
The place of lepton interaction in IC scattering, that is, the production of a photon at , after being injected at a distance from the star, at an angle (we discuss further details of geometry in the next Section and in the Appendix, especially see Fig. 23), is calculated from an inverse method. For an specific random number in the range , the interaction place is given by the formula:
(8) 
where is the rate of electron interaction to IC process (see Eq. 95 of the Appendix). From this it follows: where the integration is over the propagation path of the electron, . Knowing , we keep integrating forward until we find numerical equality of the integral result with the randomgenerated quantity , what defines the position of interaction . For a large number of simulations, the distribution of interaction distances is presented in Fig. (18, left panel).
Once we get the place of electron interaction, , simulated in the way described, we get the energy of the scattered photon from an inverse cumulative distribution. The energy of the photon produced in IC is then given by the spectrum of photons after scattering (see Appendix, Eq. 85). For a random number we define:
(9) 
where is the maximal scattered photon energy (see Appendix Eq. 72), while the energy is the MonteCarlo result for the simulated photon energy. Note that the denominator of the latter expression is the normalization needed for the MonteCarlo association to succeed. The statistics for the energies of upscattered photons is shown in Fig. 18, right panel. All relevant formulae and more details of implementation are given in the Appendix.
In a similar way we randomize by MonteCarlo the needed magnitudes for absorption. The probability of interaction for a photon with energy at a given distance, say , from an injection place (i.e., provided and are known) during the propagation in an anisotropic radiation field is given by expression:
(10) 
From this it follows, thus, similarly to the process just described, for the random number we can get the specific place of interaction . The results of such MonteCarlo simulations are presented in Fig. 19, left panel.
The energy of the lepton () produced in the photon absorption process (with the photon having an energy ) is calculated again from an inverse function. The latter is obtained from the integration of the spectra produced in the process: (see Appendix, Eq. 67):
(11) 
where is a new random number. The integration is over the electron energy . For normalizing, as the produced lepton spectra are symmetric with respect to the energy , the lower limit of this integral is fixed to this latter energy, while the upper limit is equal to maximum energy available in this process . From the energy of the electron, , the energy of the associated positron is also obtained as . Again, note that the denominator of the latter expression is the normalization needed for the MonteCarlo association to succeed.
Checks for the MonteCarlo distributions
Here we discuss the rightness of the simulated random distribution of electrons and photons resulting from the cascading process. In Fig. 20 we show the comparison between the event statistics, i.e., with being the total number of simulations run in these examples, and the analytical computed probability of interaction , where is, correspondingly, the one corresponding to absorption and ICS. We see total agreement of the MonteCarlo and analytical probabilities, i.e., whereas the position of interaction of a single photon or electron, from which the subsequent cascading process is followed, is obtained through MonteCarlo and thus is random, the overall distribution maintains the shape provided by the physical scenario: given the target and injection energy, the mean free path defines the distribution for a sufficiently high number of runs.
Parameter  Symbol  Adopted value 

Radius of star  
Mass of star  
Temperature of star  K  
Mass loss rate of star  
Wind termination velocity  
Wind initial velocity  
Distance to the system  
Eccentricity of the orbit  
Semimajor axis  
Longitude of periastron 
Basic geometry
Parameters such as the viewing angle towards the observer, , the separation of the binary, , and the distance to the shock region along the orbit, (Eq. 1), are directly connected to ray observational results and are thus crucial for detailed model discussion. Table 4 gives the set of LS 5039 orbital and binary parameters that are relevant for modeling. These parameters come from the recent work by Casares et al. (2005b), either directly from their measurements or from the values compiled by them. The inferior conjunction (INFC) is defined as the phase when the binary is viewed from the pulsar side (the smallest ). Exactly in the opposite side of the orbit we find superior conjunction (SUPC), when the pulsar is behind the massive star (the largest ), see also Fig. 1. One can see that for LS 5039, INFC () is close to apastron phase, while SUPC () is close to periastron.
Figs. 21 and 22 shows the angle to the observer, , as a function of orbital phase for two considered inclination of the binary, and . In addition, the separation of the binary is marked. The viewing angle changes within the limits and depends also on the longitude of the periastron .
The comparison for the two inclinations shows already that varies significantly in the case of the larger angle, as the difference is , what is crucial for all angledependent processes discussed in this paper. The separation of the binary is given by
(12) 
where is the eccentricity of the orbit,
(13) 
the value of is defined as
(14) 
and where is the semiminor axis. As the semimajor axis of LS 5039 is AU the assumed pulsar in the orbit is at periastron only from the massive star, whereas at apastron, it lies at , about a factor of 2 farther.
To get the relation for the angle to the observer we use the geometry shown in Fig. 21 where starting from the spherical triangle containing the inclination angle , we have the relations:
(15) 
and
(16) 
which is the cosine theorem for a rectangular spherical triangle. Then, the angle to the observer defined in Fig. 21, is calculated from the formula:
(17) 
from which:
(18) 
The sign in this relation depends on the orientation to the observer, i.e., if the pulsar is below or above the orbital plane. Note that Eq. (17) is not fulfilled all around the orbit as in here the angle denotes the common point for the orbital and observer plane. The relation can then be used after defining the orientation of the orbital plane with respect to the observer given by the inclination and the periastron position . The latter is done when we tilt the inclined orbital plane with . To find the internal phase dependent on the true anomaly (with at periastron), we have to find the spherical triangle from Fig. 21 in the system of the orbital plane. Including the tilting of the systems we have the relation:
(19) 
Then we have to check if the calculated value of is within one orbital phase . If not the angle have to be replaced by (if ) or (). To use the Eq. (17) and calculate the corresponding angle to the observer we have to recalculate the angle once again as it is only valid in the range . We can divide the orbit in four quarters in which the following transformations are needed:

if ,

if ,

if ,

if .
Inserting such phase into the Eqs. (17) and (18) we get the angle to observer corresponding to the true anomaly .
Thermal radiation
The source of the thermal radiation field is the massive star of early type (O, Be, WR). The spectrum is described by Planck’s law, which differential energy spectrum (the number of photons of given energy per unit energy , per unit solid angle , per unit volume ) is given by:
(20) 
where is thermal photon energy, is the Planck constant, and is Boltzmann constant.
Gamma ray absorption and production: opacity and geometry
The optical depth to photon absorption in the radiation field of the massive star up to infinity can then be calculated from the integral:
(21) 
where is a photon interaction rate to production in an anisotropic radiation field and is its propagation length. When the propagation occur toward the massive star surface the integration is performed up to the stellar surface. The photon interaction rate, , is related to a photon of energy injected at a distance from the massive star, at angle (see Fig. 23), and is given by the formula:
(22) 
where is the distance to the interacting photon from the injection place along its propagating path. The integration limits will be discussed in a different part of the Appendix. For now we focus on the variable in the first integration. It is the cosine of the photonphoton scattering angle (see Fig. 23). The angle is the azimuthal angle between the photon (ray) propagation direction and the direction to the massive star, while gives its limit value for the direction tangent to the massive star surface (see Fig. 25). The cross section to production is denoted as , where the parameter in the center of mass system is
(23) 
with
(24) 
being the photon energy squared (in this notation it is assumed also that ). From this we get
(25) 
The kinematic condition for the angle which defines the threshold for the creation is given by expression:
(26) 
To simplify the equations we rewrite the internal integration in Eq. (22) making use of , such that,
(27) 
With the replacement
(28) 
where
(29) 
and defining the constant , the spectrum of thermal photons is now given by the formula:
(30) 
Substituting in the internal integral , and introducing the integration variable to , via
(31) 
yields to the integral:
(32) 
The lower limit of integration is from the energy condition for the process , i.e., it follows from the threshold condition , where . . The upper integration limit comes from the relativistic limit , where we get . To proceed forward, we introduce a dumb variable, , by , to get:
(33) 
The cross section for pair production is (Jauch and Rohrlich, 1980):
(34) 
where is the classical electron radius, and is the Thomson cross section. When putting this expression into the integral (Eq. 33) we get finally,
(35) 
with .
Weparametrizee the internal integral and write it as a function of ,
(36) 
The plot for this function is shown in Fig. (24). Then, the integral we are after can be written as and
(37) 
In the second internal integral, we have to fix the limits (having in mind that the angle depends on the angle , then also on the variable ),
(38) 
The angle determines the maximal azimuthal angle of ray photon propagation with respect to the direction of the thermal photon, so that it gives the directions tangent to the star surface. This condition for can be determined from the spherical triangle (shown in Fig. 25). From the cosine theorem of spherical trigonometry
(39) 
where is defined as in Fig. (25). After some algebra we get:
(40) 
where , while . The angle determines the direction of incoming ray photon (see Fig. 23).
The limits of integration with respect to the parameter (see Fig. 25), are the range of angles for soft photons coming from the star,
(41) 
where . Using basic trigonometrical dependencies between the angles, and (Fig. 25, left), we get two expressions:
(42) 