Probing orbital parameters of gamma-ray binaries with TeV light curves

Gamma-ray binaries are binary systems where the energy flux peaks in the gamma-ray energy band. They harbour a compact object (a neutron star or a black hole) orbiting around a massive star that provides a strong radiation field. It is believed that the gamma-ray emission from such objects can be strongly attenuated through the electron-positron pair production in gamma-gamma interactions. The importance of gamma-gamma absorption depends on the orbital phase and on the geometry of the system. In this work we propose a method of how the orbital parameters of gamma-ray binaries could be probed with TeV light curves that have imprinted features of gamma-gamma absorption.


INTRODUCTION
Gamma-ray binaries constitute a small but growing class of high mass binary systems with spectral energy distributions which peak above 1 MeV in a νF ν representation (see Dubus 2013Dubus , 2015;;Chernyakova et al. 2019;Chernyakova & Malyshev 2020, for a detailed review of these systems).The class comprises only nine known objects: eight in our Galaxy and one in the Large Magellanic Cloud.These systems consist of a compact source (pulsar or black hole) orbiting around a massive companion star and for only two of them is the nature of the compact source well established -both PSR B1259-63/LS 2883 and PSR J2032+4127/MT91 213 comprise of a pulsar.More recently, radio pulsations have been reported for LS I +61 303 Weng et al. (2022).While PSR B1259-63/LS 2883 is a very well studied binary, PSR J2032+4127/MT91 213 was confirmed as a gamma binary more recently, due to its very long orbital period of 45 − 50 years.Very-high-energy gamma-ray emission was detected in 2017 as it passed periastron (Abeysekara et al. 2018).
Gamma-ray emission in such plerionic binaries is believed to be produced through the inverse Compton scattering of iurii.sushch@nwu.ac.za vansoelenb@ufs.ac.za electrons, accelerated at the termination shock that forms between the pulsar and stellar winds, on the stellar radiation field and thus is strongly orbital dependant.Additionally, hydrodynamic simulations imply the formation of the Coriolis 1 shock and other secondary shocks (see e.g.Bosch-Ramon et al. 2015;Dubus et al. 2015;Huber et al. 2021a) resulting in a full confinement of particles and providing additional acceleration and re-acceleration sites.Simulations suggest the highest density of very-high-energy electrons and hence the strongest gamma-ray emission at the apex of the shock (Dubus et al. 2015;Huber et al. 2021a) but this result relies on certain assumptions, e.g.prescription for the particle injection which could be oversimplified.Therefore, although the most plausible location for the emission region is the apex of the shock, emission can be produced at the other sites, located also behind the pulsar with respect to the star, and it is not trivial to address the significance of their contribution.
For the accretion powered systems with a black hole as a compact object (known as microquasars) particles can be accelerated anywhere along the jet through various mechanisms (see e.g.Dubus 2006a; Khangulyan et al. 2008;Bosch-Ramon & Khangulyan 2009).Recent H.E.S.S. observations of the SS 433 system provide for the first time a direct mea-1 Note, that Dubus et al. (2015) argue that the existence of the "back shock" is not related to Coriolis forces as it also forms in the absence of rotation arXiv:2310.01102v1[astro-ph.HE] 2 Oct 2023 surement of the location of the particle acceleration site in a jet-like source which appears to be at a considerable distance of ∼ 25 pc from the central engine (Olivera Nieto 2023).Therefore, the gamma-ray emission in these systems can be produced far off the orbital plane and far from the compact object and its companion star.Gamma-ray photons emitted in binary systems are subject to γγ absorption as they travel through the photon field created by the massive star (see e.g.Dubus 2006b).Moreover, γγ absorption might be the main reason for the characteristic minimum in the TeV light curve observed in gamma-ray binaries (Dubus 2006b;Sushch & van Soelen 2017).In the case of PSR B1259-63/LS 2883, a very eccentric binary (e = 0.87) with an orbital period of 3.4 years, the TeV flux from the system increases as the pulsar moves closer to the star but then suddenly drops as it gets closer to periastron itself and then increases again after the periastron passage forming the second peak 2 (Aharonian et al. 2005(Aharonian et al. , 2009;;H.E.S.S. Collaboration et al. 2013;Romoli et al. 2015;H. E. S. S. Collaboration et al. 2020).This behavior is counter-intuitive, because at periastron the pulsar is located at the closest distance to the star and thus encounters the highest density of the target photon field for inverse Compton radiation.However, these conditions are also optimal for gamma-gamma absorption which might cause a severe decrease of the TeV flux.It was shown for PSR B1259-63/LS 2883 that the orbital phase for which the absorption would be the strongest coincides with the orbital phase of the dip in the light curve (Sushch & van Soelen 2017) hinting that this decrease of the TeV flux indeed could be related to the gamma-gamma absorption.The level of absorption is, however, not sufficient to explain the magnitude of the dip in the TeV light curve if the emitting region is located close to the pulsar.Nevertheless, it is interesting to note that in the phase-folded stacking analysis presented in H. E. S. S. Collaboration et al. (2020) there is a hint of a hardening of the TeV spectrum before the periastron (∼ −20 d) which is coincident with a dip in the light curve.This is what would be expected from a change induced by significant gamma-gamma absorption as proposed in Sushch & van Soelen (2017).For LS 5039, a much more compact binary with low eccentricity and an orbital period of only 3.9 days, the modulation of the light curve could be well explained by gamma-gamma absorption with the opacity at the dip high enough to explain the level of flux attenuation (Dubus 2006b;Böttcher & Dermer 2005;Huber et al. 2021b).

2
It should be noted that there is some evidence for a less prominent third "middle" peak right before the periastron passage (see the analysis of the H.E.S.S. data for all periastra passages since 2004 in H. E. S. S. Collaboration et al. 2020) If the emitting region is located close enough to the star, i.e. the gamma-gamma opacity is high enough to modulate the TeV light curve, the location of the minimum in the light curve would contain information on the geometry of the system, because such orbital parameters as the inclination angle and the longitude of periastron determine the orbital phase where the absorption would be the highest.
In this paper we show how, if it is assumed that gammagamma absorption is significant enough to attenuate the very high energy emission, the location of the minimum in the light curve can be used to place constraints on the orbital geometries the gamma-ray binaries.This is applied to the system PSR J2032+4127/MT91 213.

Gamma-gamma absorption
The interaction of a γ-ray photon of energy ϵ γ with a lowenergy photon of energy ϵ, can result in electron-positron pair production (Gould & Schréder 1967), if the energy exceeds the threshold condition, given by where the energies are normalized to the electron rest-mass energy, i.e. ϵ = hν/m e c 2 , and θ int is the interaction angle between the two photons.
The γγ optical depth is given by (Gould & Schréder 1967) where l is the distance over which the γ-ray photon travels, µ = cos θ int , dΩ = sin θdθdϕ is the solid angle element in the spherical coordinate system centered at the gamma-ray photon with zenith determined by the direction from the star (see Appendix B), and n ph (ϵ, Ω) is the number density per unit solid angle of the low-energy target photons.Here, σ γγ , is the γγ cross-section (Jauch & Rohrlich 1976), where and σ T is the Thomson cross-section.The radiation field, which provides the target photons for inverse Compton scattering and γγ absorption, is produced by the massive companion star.In the case of a Be star companion, the photon field is produced by two components: the optical stellar radiation and the infrared radiation from the circumstellar disc.Calculations of the γγ absorption in PSR B1259-63/LS 2883 due to the circumstellar disc (Sushch & van Soelen 2017) showed that the circumstellar disc component is considerably less important than the stellar radiation field.The energy densities of the two radiation fields are constrained by the observed luminosities of the star and of the circumstellar disc.Although the energy density of the circumstellar disc photons can be higher the effective interaction length is much smaller and hence the contribution to the overall γγ absorption is insignificant.In this study for the simplification of the method we only take into account stellar photons.We then approximate the photon density distribution by black-body radiation and assume that photons are emitted radially.Therefore, where h is the Planck constant, c is the speed of light, T * is the effective temperature of the star, and B ν (T * ) is the Planck function, i.e. the spectral power emitted per unit area per unit solid angle.
It is clear from Eq. 2 that for a known energy and angular distribution of target photons the optical depth would depend exclusively on the energy and trajectory of a gamma-ray photon.The trajectory of a gamma-ray photon traveling along the line of sight in turn depends on the the geometry of the system and the location of the emitting region at a particular moment of time.The location of the emitting region or region characterized by the most effective particle acceleration in gamma-ray binaries is unclear and is still under active debate.The emitting location will be different depending on if the compact object is a pulsar or a black hole.In the former case particles are believed to be accelerated at the termination shock between the pulsar and stellar winds as well as at other shocks generated due to the winds interaction (see e.g.Dubus et al. 2015;H. E. S. S. Collaboration et al. 2020;Huber et al. 2021a), while in the latter case particles can be accelerated anywhere along the jet (see e.g.Dermer & Böttcher 2006;Bosch-Ramon & Khangulyan 2009).Moreover, it is not completely clear which part of the termination shock (or other shock) in the wind-driven scenario would be the most favourable for particle acceleration and this can also change with orbital phase of the pulsar (see Section 1).It should be noted that this can have a dramatic impact on the expected emission spectrum as well as the expected gamma-gamma absorption.In the wind-driven scenario the apex of the bow shock is argued to dominate the gamma-ray emission.In this case γγ absorption could have a strong effect due to the proximity to the companion star.But if the gamma-ray emission is produced at the back (Coriolis) or secondary shocks the emission region could be too far from the star for any significant pair production to take place.Similarly, for microquasars, the effect of γγ absorption would be very different depending on whether the emission is produced at the base of the jet or further along the jet.
In the following we eliminate this effect from our considerations by assuming a point-like scenario, i.e. that particles are accelerated, and hence the gamma-ray emission is produced, close to the compact object.We acknowledge that this assumption would be realistic only for the accretion powered systems where gamma-ray emission is generated at the base of the jet being launched close to the central engine.For the purposes of the proposed method this assumption is, however, also well justified for wind driven systems with the pulsar as a compact object.Depending on where exactly the dominant fraction of the gamma-ray emission is produced (apex or tail of the shock) the level of the γγ absorption will change, but it will not significantly impact the orbital dependence of the optical depth and specifically the orbital phase of the maximum absorption, which is of prime interest for the proposed method, because the location of the emitting region is in (or close to) the orbital plane.This was demonstrated for the case of PSR B1259-63/LS 2883 in Sushch & van Soelen (2017) where the optical depth was calculated both assuming the gamma-ray emission is produced at the location of the pulsar and at the location of the apex of the shock.In some systems the massive companion star features a circumstellar disc which would push the wind termination shock closer to the pulsar during disc crossings.In this case the apex of the shock would actually be closer to the location of the pulsar.The situation can be more complicated if the dominant site of the gamma-ray emission production is changing with the orbital phase.In this case the proposed method would not be applicable as the dependence of the optical depth on the orbital phase can be regulated by how the location of the emission region is changing.This scenario is, however, not well understood and beyond the scope of this paper.Finally, this method would also not be applicable for the microquasar scenario where the gamma-ray emission is produced along the jet at a considerable distance from the black hole.In this case the dependence of the γγ absorption on the orbital phase can be very different and would depend on the orientation of the jet.
Adopting the assumptions above, we end up with the optical depth being exclusively a function of the geometry of the system, i.e. the inclination angle i (the angle between line-of-sight (LoS) and the normal to the orbital plane), the longitude of periastron ω (the angle between the ascending node and periastron), the eccentricity e, and the angular orbital phase χ which we choose to measure from apastron in the direction of motion of the compact object (true anomaly minus π; see Fig. 1).Further, following the hypothesis that the orbital phase of the characteristic dip in the TeV light curve of a binary is defined by the strongest absorption, we can probe the geometry of the system. .Schematic cartoon depicting the geometry of the binary system projected on the orbital plane.The LoS direction shows the line-ofsight towards the observer, χ is the angular orbital phase measured from apastron and ω is the longitude of periastron.For clarity we also show on the plot what the value of ω would be if the LoS direction is aligned with the blue arrows.

Analytic considerations
Before performing full numerical calculations, we present the derivation of a simple analytic approximation of the γγ optical depth as a function of the orbital phase χ measured from apastron in the direction of the movement (see Fig. 1).For this we approximate the γγ cross section by (Zdziarski & Lightman 1985): eliminating the dependence of the cross section on the energy of the soft photon and the interaction angle.Figure 3.12 in Boettcher et al. (2012) shows the comparison of the pair production rate spectra using different approximations for the cross section.The delta-function approximation works very well for intermediate energies of resulting pairs but predicts artificially hard cutoffs compared to the full expression of the cross section.Assuming the stellar black body radiation as the main target photon field and adopting a point-like approximation we can express the number density of the photon field as where r is the distance from the gamma-ray photon to the star and x is the distance the gamma-ray photon has moved along the line of sight.The delta-function reflects the assumption that the stellar photons are radiated radially from the star taking into account that the zenith is aligned with the direction from the star to the gamma-ray photon as described above.
Then, taking into account both delta-functions Eq. 2 can be rewritten as where is only a function of the gamma-ray energy and is constant with the orbital phase.The function reflects the geometry of the system and determines the variability of the optical depth with the orbital phase.The terms µ(x, χ) and r(x, χ) can be expressed as functions of the distance between the star and gamma-ray emitting region d and the initial interaction angle, µ 0 ≡ µ(x = 0), as As long as the observer is far enough away from the binary system, the integration in Eq. 10 can be performed up to infinity.Then substituting Eqs.11 and 12 f (χ) in Eq. 10 and performing the integration we find that f (χ) can be expressed analytically as (see Appendix A for details) The maximum of this function determines the orbital phase where the maximum gamma-gamma absorption occurs and hence, according to our hypothesis, the orbital phase of the dip in the TeV light curve.Now, applying a point-like source approximation, i.e. the emission region is located at the position of the compact object, we can express d(χ) and µ 0 (χ) in terms of the orbital parameters of the system: Following these simple considerations it can be seen that the maximum of the function f (χ, i, ω, e) in Eq. 13 is determined by the geometry of the binary system, namely its inclination angle, longitude of periastron, and eccentricity.Below we will demonstrate that this approximation is in fairly good agreement with the exact numerical calculations for a wide parameter space and hence can be used to place rough constraints on the geometry of the binary system.

Comparison of analytic and numeric solutions
Numeric calculation of the optical depth were done following the approach presented in (Sushch & van Soelen 2017) which for the stellar radiation component essentially follows Dubus (2006b).The method is summarized in the Appendix B. To explore the deviation of the analytic approximation from the exact numeric simulation we consider a Test Binary with parameters listed in the Table 1.The optical depth is calculated both analytically and numerically varying separately the inclination, the longitude of periastron and the eccentricity.Other parameters are kept fixed.In Figure 2 we show the γγ optical depth τ γγ as a function of the orbital phase χ with color coding indicating the range of varied parameters: left panels -varied inclination with the longitude of periastron fixed at 135 • and eccentricity fixed at 0.5; middle panels -varied longitude of periastron with the inclination fixed at 45 • and the eccentricity fixed at 0.5; and right panels -varied eccentricity with the inclination fixed at 45 • and the longitude of periastron fixed at 135 • .The upper panels show the analytical approximation discussed above, the middle panels show the exact numeric solutions of Eq. 2, and the lower panels show direct comparisons for some selected parameters.
In general, there is a rather good agreement between the analytic and numeric solutions.Both solutions show the same evolution of the shape of the optical depth curve with the change of the orbital parameters.The analytic approximation slightly overestimates the value of the optical depth which becomes more significant at large inclination angles.This, however, does not play a major role for the purpose of this paper.
A careful parameter scan showed that the discrepancy between the analytic approximation and numeric calculations in determination of the orbital phase of maximum absorption becomes more significant for a specific case of ω ∼ 270 • , i.e. the periastron coincides with the inferior conjunction.This is particularly true for high eccentricities where the optical depth curve exhibits two maxima for ω ∼ 270 • (Fig. 3).
Figure 4 shows the application of the optical depth calculation to two well studied binary systems, PSR B1259−63 (top panel) and LS 5039 (bottom panel), with substantially different orbital parameters (Table 1).While PSR B1259 − 63 is a very eccentric and long period binary, LS 5039 is a compact binary with a low eccentricity.In both cases the analytic approximation reproduces the shape of the numerically calculated curve relatively well with a good match of the orbital phase of the peak optical depth.

The orbital phase of the maximum absorption
To estimate the orbital phase where the gamma-gamma absorption is the strongest and hence where the dip in the observed light curve is expected we vary the geometrical parameters of the orbit, namely inclination between 0 and 90 degrees (20 bins) and longitude of periastron between 0 and 360 degrees (20 bins), and for each combination of i and ω we calculate the orbital phase with the highest optical depth, i.e. maximum absorption, ϕ max .For the results presented below we use exact numeric simulations but the analytical approximation is compatible with these results.The dependence of the maximum absorption phase on the orbital parameters is further discussed and a detailed comparison of the numeric and analytic approach is presented in Appendix C. For numeric calculation we split the orbit into 200 bins equidistant in orbital phase setting a systematic error of of 1.8 • on the estimate of ϕ max .
3. APPLICATION TO PSR J2032+4127/MT91 213 PSR J2032 + 4127/MT91 213 is only the second known gamma-ray binary for which it is well established that the compact object is a pulsar (Abeysekara et al. 2018).The system consists of the pulsar PSR J2032 + 4127, which was first detected at GeV energies by Fermi-LAT (Abdo et al. 2009) and later confirmed in radio observations (Camilo et al. 2009), and a ∼ 15M ⊙ Be star MT91213 (Lyne et al. 2015).The periastron passage in 2017 presented us with a unique opportunity to detect this binary system at VHE which otherwise would be impossible given its extremely eccentric orbit   VHE observations of the PSR J2032 + 4127/MT91 213 binary system around the 2017 periastron passage were conducted by both VERITAS and MAGIC (Abeysekara et al. 2018) and revealed a firm detection of the variable emission associated with the binary.The observed TeV light curve shows a rather typical behaviour with a two-bump structure and a dip shortly after periastron.It is interesting that the shape of the X-ray light curve is very different from the TeV light curve.It exhibits a much broader dip with the dimming starting already ∼ 40 days before periastron and the minimum flux period spanning from shortly before to 5 − 10 days after periastron (Abeysekara et al. 2018).Although the TeV dip is coincident with the minimal X-ray flux, so is the first TeV bump which occurs roughly at periastron where the Xray flux is already at its minimum.These differences indicate that the decrease of X-ray and gamma-ray flux could be attributed to different processes which agrees with the gammagamma absorption scenario for the attenuation of the TeV flux.TeV gamma-ray photon as a function of the orbital phase from the binary system PSR B1259 − 63 (top panel) and LS 5039 (bottom panel; orbital solution by Casares et al. 2005).Filled orange circles represent the exact numeric estimate of Eq. 2, while the solid blue line shows the analytic approximation derived in Section 2.2.
To probe the orbital parameters of PSR J2032+4127/MT91 213 we follow the procedure described in Section 2.4.In Fig. 5 we show numerically calculated maps depicting the orbital phase where the maximum absorption occurs (top panel) and the value of the optical depth at that orbital phase (bottom) for the full parameter space of i and ω.The spread of potential maximum absorption locations covers almost the whole orbit depending on the inclination and the longitude of periastron.However, only for high inclincation angles is the optical depth high enough to cause significant absorption.It should be noted here that the assumption of the emitting region being located at the pulsar position might somewhat underestimate the optical depth.For these simulations we used an eccentricity of e = 0.961 and a period of P = 17000 days, which is the average 'model 2' presented in Ho et al. (2017).According to simulations presented in Fig. 2 small variations in the eccentricity do not strongly impact the estimate of the orbital phase with the maximum absorption.They can, however, significantly change the estimate of the optical depth itself, resulting in a stronger absorption at ϕ max for higher eccentricity.The estimate of ϕ max is not very sensitive to the orbital period as long as it is long enough and the size of the star is negligible (see Appendix C).For other parameters fixed, a short orbital period would result in stronger absorption as the distance between the pulsar and the star is smaller.
Both MAGIC and VERITAS observations (Abeysekara et al. 2018) indicate that the minimum in the light curve occurs at about 7 days after the periastron passage.In our diagnostic we conservatively consider a larger time interval of 8 days spanning from 3 days to 11 days after periastron.Conversion of days to the orbital phase is naturally very sensitive to the eccentricity and the orbital period, estimates of which suffer from quite large uncertainties, for PSR J2032 + 4127/MT91 213.
Taking into account the whole range of values of e and P allowed by optical observations (see Table 1; Ho et al. 2017) we end up with quite a large range of orbital phases, 186 • − 290 • , which could correspond to the minimum in the light curve.Nevertheless, even this rather wide range of orbital phases is still quite constraining in the determination of the inclination angle and the longitude of periastron.In Fig. 6 the regions marked with green colors represent the parameter space allowed if the dip in the TeV emission is due to gammagamma absorption, i.e. those that fall into the 186 • − 290 • range, as implied by numeric calculations (top panel) and analytic approximation (bottom panel).One can see that it covers roughly half of the whole parameter space.However, taking into account that the optical depth needs to be high enough to provide a sufficient level of absorption, one can further strongly reduce the parameter space.The dark green color corresponds to the parameter space where at least half to the intrinsic radiation is absorbed.It is clear that taking into account the level of absorption strongly constrains the inclination, suggesting a highly inclined orbit.Naturally a better estimate of the orbital period and eccentricity further constrains the allowed parameter space.The orange region shows the parameters allowed by 'model 2' of (Ho et al. 2017) that adopts average values of P = 17000 days and e = 0.961.This allows us to much better constrain the longitude of periastron.Again, the dark orange region reflects significant levels of absorption.
The dashed regions in Fig. 6 illustrate the constraints on i and ω as set by optical observations (Ho et al. 2017).Their comparison to the constraints obtained in this work can be summarized in a few main points: • a clear overlap of the regions with allowed values strongly supports the general idea that the location of the minimum in the TeV light curve is determined by the gamma-gamma absorption and reflects the geometry of the system • a precise measurement of the orbital phase of the dip in the TeV light curve, that also requires a good knowledge of the orbital period and eccentricity, could potentially not only strongly constrain the i−ω parameter space but also yield quite precise estimates of at least one of these two parameters, namely the longitude of periastron.
• the estimate of the optical depth which is a measure of the level of absorption might strongly constrain the inclination of the system under condition of good understanding of where the emitting region is located.
• vice versa, in the case of reliable estimates of the orbital geometry from other considerations like optical observations or radio/gamma-ray pulsar timing measurements, the location of the dip in the TeV light curve can indirectly indicate the location of the emitting region and hence provide a valuable insight into understanding of the particle acceleration in these systems.
The latter is certainly the case for the PSR B1259-63/LS 2883 binary and will be further investigated in our future works.
Comparison of the top (numeric calculations) and bottom (analytic approximation) panels in Fig. 6 does not show a big difference between the two approaches implying that the analytic approximation could be a useful tool that can be used to constrain the geometry of the orbit without spending a lot of CPU hours.

SUMMARY
In this study we examine the impact of gamma-gamma absorption on the periodic very-high-energy radiation from gamma-ray binaries and study how this impact depends on the orbital parameters of a particular system.We suggest that under the assumption that a dip-like feature in the TeV light curve (characterized by a fast decrease flux superposed on a general increase of flux) can be attributed to the highest gamma-gamma absorption one can infer constraints on orbital parameters of the system from the orbital location of this dip.We propose a method to diagnose TeV light curves and use them as a probe for orbital parameters.We also derive an analytic approximated solution for the optical depth that agrees reasonably well with the exact numeric solution and provides an efficient tool for quick diagnostics.Application to PSR J2032 + 4127/MT91 213 results in constraints on the inclination and longitude of periastron that are in a good agreement with constraints obtained from optical observations, implying that indeed this method could be used for characterization of the orbit.Moreover, we argue that precise time-resolved flux measurements as well as accounting for the level of absorption would further constrain the orbital parameters.Therefore, the method of gamma-gamma absorption offers another completely independent way of probing orbital parameters that could be complementary to classic optical radial velocity measurements.On the other hand, well determined orbital geometry combined with the TeV light curve could point to the location of the emitting region and shed light on acceleration and emission processes in these systems.Substituting Eqs.11 and 12 f (ϕ) in Eq. 8 and changing the upper limit of integration to infinity we get We further substitute y = x + dµ 0 and split it into two integrals For the first integral I 1 , .

B. NUMERIC SIMULATIONS
The integration over a solid angle in Eq. 2 is performed in the spherical coordinate system centered at the location of the emitting region which is assumed to overlap with the location of the compact object.The zenith is determined by the direction from the center of the star to the emitting region.The integration in zenith angle can be then substituted by the integration in η = cos θ and limited to η ∈ [cos θ * ; 1], with θ * being the angular radius of the star as viewed from the gamma-ray location, i.e.
, where d is the distance between the gamma-ray location and the center of the star and R * is the radius of the star.Additionally, we set the azimuth angle to be measured from the projected direction towards the observer and integrate over ϕ from 0 to π using the symmetry.The solid angle integral in Eq. 2 can be then rewritten as where η * = cos θ * .The integral along the line-of-sight is calculated in cylindrical coordinates centered on the center of the star and tied to the orbital plane from the location of the emitting region up to 1000 stellar radii.Note, that in Sushch & van Soelen (2017) the spherical coordinate system in which the integration over the solid angle is performed is defined differently due to accounting for the additional radiation field supplied by the circumstellar disk.It is centered at the emitting region but with zenith aligned with the normal to the disk.In calculations performed in Sushch & van Soelen (2017) we mistakenly flipped the direction of the soft photon to the opposite direction, i.e. the photon radiated from center of the star was calculated as the photon traveling towards the star, which resulted in miscalculation of the interaction angle at every step of the integral along the line-of-sight and eventually in overestimation of the optical depth both for the stellar radiation and the circumstellar disc (Fig. 7).This mistake does not considerably change the shape of the orbital evolution of the optical depth.Qualitatively the results presented in Sushch & van Soelen (2017) are not strongly affected, but the overall effect of the gamma-gamma absorption is expected to be weaker.All the results will be updated in the forthcoming dedicated paper on PSR B1259-63/LS 2883.

C. COMPARISON OF ANALYTIC AND NUMERIC APPROACH
To calculate the orbital phase with the maximum absorption, ϕ max , we split the orbit into 200 bins setting a systematic error of 1.8 • .We estimate ϕ max for different combinations of i and ω varying the inclination between 0 and 90 degrees (20 bins) and longitude of periastron between 0 and 360 degrees (20 bins).Figures 8 and 9 show maps depicting color-coded ϕ max as a function of i (y-axis) and ω (x-axis) calculated numerically (top panels) and analytically (upper middle panels).The lower middle panels show the difference between two estimates, ϕ num max − ϕ an max , and the bottom panels map the the maximum optical depth as obtained from numerical simulations, i.e. the optical depth at ϕ num max .The stellar parameters used for these simulations correspond to the Test Binary entry in the Table 1. Figure 8 explores the dependence on the eccentricity showing results for e = 0.1 (left column), e = 0.45 (middle column), and e = 0.9 (right column) for the same orbital period of 100 year, while Figure 9 shows the dependence on the period for the same eccentricity of e = 0.5 with P = 10 days (left column), P = 100 days (middle column), and P = 10, 000 days (right column).
As expected ϕ max strongly depends on the longitude of periastron of the system as long as the system is sufficiently inclined, while for low inclination angles (observer looks at the system face-on) ω is irrelevant for the estimate of ϕ max as it stays at ≈ 180 • (periastron) determined by the to the smallest binary separation at this phase.For i ≳ 10 • , ϕ max does not strongly depend on i for the orbits with low eccentricity and is mainly determined by the longitude of periastron For high eccentricities the dependence on the inclination angle becomes more significant in the domain where ω is close to 270 • (periastron coincides with inferior conjunction).This is an exclusively geometrical effect reflecting the fact that for high eccentricity orbits the emitted gamma-ray photon would be traveling closer to the star when emitted at some offset to the apastron and hence the maximum absorption would occur at some orbital phase either before or after the periastron.Additionally, as mentioned previously for high eccentricity orbits the optical depth curve for ω ∼ 270 • features two maxima that further complicates things.
Apart for a small domain around ω ≃ 270 • and the case with a very low orbital period, numeric and analytic results are in good agreement with the difference at the level of the systematical error of the numerical calculations.In the domain around ω ≃ 270 • the analytic solution becomes significantly different (although still at most by ∼ 15 • ) and it becomes more relevant at high eccentricities as this uncertainty region covers a larger range of inclination angles.However, this part of the parameter space also corresponds to weaker absorption as can be seen from the bottom panels which maps τ γγ (ϕ num max ).Indeed, in this case the difference between the numeric and analytic solution can be irrelevant as the absorption might simply be not strong enough to modify the observed light curve and hence the proposed method cannot be applied.Note, that for the high eccentricity case (right column in Fig. 8) although the the difference between the numeric and analytic solutions does correlate with the weakness of the absorption, the optical depth is high for the whole parameter space.This is due to the chosen parameters as for the e = 0.9 case the orbital period of 100 days corresponds to a very small separation distance between the pulsar and the star at periastron.These compact binary systems, however, are also irrelevant to the proposed method simply because most of the gamma-ray emission will be absorbed at any orbital phase and therefore not detected.
The numeric estimate of ϕ max is not very sensitive to the orbital period as long as it is not too short.There is no obvious difference between the maps for P = 100 days and P = 10000 days, which is confirmed by almost identical "difference" maps.Note, all three maps depicting the analytic solutions are identical as the analytical solution is independent of the orbital period (see Section 2.2).For very short periods, i.e. compact, systems (see P = 10 days column) the difference becomes more apparent, mainly due to a non-negligible size of the star.Nevertheless, as already mentioned above, these very compact systems are not of interest to us due to the impossibility of their detection.
Figure1.Schematic cartoon depicting the geometry of the binary system projected on the orbital plane.The LoS direction shows the line-ofsight towards the observer, χ is the angular orbital phase measured from apastron and ω is the longitude of periastron.For clarity we also show on the plot what the value of ω would be if the LoS direction is aligned with the blue arrows.

Figure 2 .
Figure 2. Optical depth τ γγ as a function of the orbital phase calculated analytically (top panels) and numerically (middle panels) for varying values of inclination (left panels) longitude of periastron (middle panels) and eccentricity (right panels).A direct comparison of the optical depth calculated analytically (dashed lines) and numerically (solid lines) for different values of corresponding parameters is shown in the bottom panels.

Figure 3 .
Figure 3. Same as the right bottom plot in Fig. 2 but for ω = 270 • .

Figure 4 .
Figure4.The optical depth for gamma-gamma absorption of the 1 TeV gamma-ray photon as a function of the orbital phase from the binary system PSR B1259 − 63 (top panel) and LS 5039 (bottom panel; orbital solution byCasares et al. 2005).Filled orange circles represent the exact numeric estimate of Eq. 2, while the solid blue line shows the analytic approximation derived in Section 2.2.

Figure 5 .
Figure 5.The color maps shows the orbital phase of the pulsar at which the emitted gamma-ray photon would encounter the highest gamma-gamma absorption along the line of sight as a function of assumed inclination and longitude of periastron (top panel) as well as the value of the maximum optical depth (bottom panel).

Figure 6 .
Figure6.Allowed combinations of values of i and ω as derived from the position of the minimum in the TeV light curve of PSR J2032 + 4127/MT91 213 that takes place at 7 ± 4 days after the periastron using numeric calculations (top panel) and analytic approximation (bottom panel).The green colors take into account the whole range of uncertainties for the estimates of the eccentricity and the period while orange colors assume the 'model 2' values of e = 0.961 and P = 17000 days from(Ho et al. 2017).Dark orange and dark green correspond to the region where the optical depth is high enough to attenuate the intrinsic flux by at least a factor of 2. Light yellow hatching represents the constraints on i and ω set by optical observationsHo et al. (2017) .

Figure 7 .
Figure 7.The optical depth for gamma-gamma absorption of the 1 TeV gamma-ray photon by the stellar radiation field (left panel) and circumstellar disc (right panel) as a function of days from periastron for the binary system PSR B1259 − 63 as calculated in Sushch & van Soelen (2017) with a geometric mistake (orange dashed line) and corrected in this work (blue filled circles).Note that the calculation of optical depth due the stellar photons accounts for the absorption of the stellar radiation by the circumstellar disc.

Figure 8 .
Figure 8.The upper two rows show the maps of the orbital phase of maximum absorption, phi max , color-coded from 0 • to 360 • , calculated numerically (top panels) and analytically (upper middle panels) for the Test Binary for three different values of eccentricity and the orbital period of 100 days.The lower middle row shows the mapped difference, ϕ num max − ϕ an max , and the bottom maps show numerically calculated maximum optical depth, τ γγ (ϕ num max ).

Figure 9 .
Figure 9. Same as Fig. 8 but for different values of the orbital period and e = 0.5

Table 1 .
Physical parameters of examined systems.