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

Iurii Sushch Centre for Space Research, North-West University, 2520
Potchefstroom, South Africa  Astronomical Observatory of Ivan Franko National
University of Lviv, Kyryla i Methodia 8, 79005 Lviv, Ukraine
<EMAIL_ADDRESS>Brian van Soelen Department of Physics, University of
the Free State, PO Box 339, Bloemfontein 9300, South Africa
<EMAIL_ADDRESS>
###### Abstract

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.

Gamma-ray sources (633) — Pulsars (1306) — Non-thermal radiation sources
(1119) — Binary stars (154) — Companion stars (291)

##  1 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 $\nu
F_{\nu}$ representation (see Dubus, 2013, 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 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 Coriolis111Note, 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 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 measurement of the location of the
particle acceleration site in a jet-like source which appears to be at a
considerable distance of $\sim 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 $\gamma\gamma$
absorption as they travel through the photon field created by the massive star
(see e.g. Dubus, 2006b). Moreover, $\gamma\gamma$ 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 peak222 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) (Aharonian et al., 2005, 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 ($\sim-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).

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 gamma-gamma 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.

##  2 Method

###  2.1 Gamma-gamma absorption

The interaction of a $\gamma$-ray photon of energy $\epsilon_{\gamma}$ with a
low-energy photon of energy $\epsilon$, can result in electron-positron pair
production (Gould & Schréder, 1967), if the energy exceeds the threshold
condition, given by

$\epsilon\epsilon_{\gamma}(1-\cos\theta_{\mathrm{int}})\geq 2,$ (1)

where the energies are normalized to the electron rest-mass energy, i.e.
$\epsilon=h\nu/m_{\rm e}c^{2}$, and $\theta_{\mathrm{int}}$ is the interaction
angle between the two photons.

The $\gamma\gamma$ optical depth is given by (Gould & Schréder, 1967)

$\tau_{\gamma\gamma}=\int_{0}^{l}\,{\rm d}l\int_{4\pi}{\rm
d}\Omega\,\,(1-\mu)\int_{\frac{2}{\epsilon_{\gamma}(1-\mu)}}^{\infty}{\rm
d}\epsilon\,n_{\rm
ph}(\epsilon,\Omega)\sigma_{\gamma\gamma}(\epsilon,\epsilon_{\gamma},\mu)$ (2)

where $l$ is the distance over which the $\gamma$-ray photon travels,
$\mu=\cos\theta_{\mathrm{int}}$, ${\rm d}\Omega=\sin{\theta}d\theta d\phi$ 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_{\rm ph}(\epsilon,\Omega)$ is the number density per unit
solid angle of the low-energy target photons. Here, $\sigma_{\gamma\gamma}$,
is the $\gamma\gamma$ cross-section (Jauch & Rohrlich, 1976),

$\sigma_{\gamma\gamma}(\beta)=\frac{3}{16}\sigma_{\rm
T}(1-\beta^{2})\left[(3-\beta^{4})\ln\left(\frac{1+\beta}{1-\beta}\right)-2\beta(2-\beta^{2})\right],$
(3)

where

$\beta=\sqrt{1-\frac{2}{\epsilon\epsilon_{\gamma}(1-\mu)}},$ (4)

and $\sigma_{\rm T}$ is the Thomson cross-section.

The radiation field, which provides the target photons for inverse Compton
scattering and $\gamma\gamma$ 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 $\gamma\gamma$ 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 $\gamma\gamma$ 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,

$n_{\rm ph}(\nu,\Omega)=\frac{1}{h\nu c}B_{\nu}(T_{\ast}),$ (5)

where $h$ is the Planck constant, $c$ is the speed of light, $T_{\ast}$ is the
effective temperature of the star, and $B_{\nu}(T_{\ast})$ 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
$\gamma\gamma$ 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 $\gamma\gamma$ 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 $\gamma\gamma$ 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 $\gamma\gamma$ 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 $\omega$ (the angle between the ascending
node and periastron), the eccentricity $e$, and the angular orbital phase
$\chi$ which we choose to measure from apastron in the direction of motion of
the compact object (true anomaly minus $\pi$; 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.

Figure 1: Schematic cartoon depicting the geometry of the binary system
projected on the orbital plane. The LoS direction shows the line-of-sight
towards the observer, $\chi$ is the angular orbital phase measured from
apastron and $\omega$ is the longitude of periastron. For clarity we also show
on the plot what the value of $\omega$ would be if the LoS direction is
aligned with the blue arrows.

###  2.2 Analytic considerations

Before performing full numerical calculations, we present the derivation of a
simple analytic approximation of the $\gamma\gamma$ optical depth as a
function of the orbital phase $\chi$ measured from apastron in the direction
of the movement (see Fig. 1). For this we approximate the $\gamma\gamma$ cross
section by (Zdziarski & Lightman, 1985):

$\sigma_{\gamma\gamma}=\frac{1}{3}\sigma_{\rm
T}\epsilon\,\delta(\epsilon-\frac{2}{\epsilon_{\gamma}}).$ (6)

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

$n_{\rm ph}\left(\epsilon,\theta\right)=m_{\rm
e}c^{2}\frac{8\pi}{h^{3}c^{3}}\left(\epsilon m_{\rm
e}c^{2}\right)^{2}\frac{1}{e^{\epsilon m_{\rm
e}c^{2}/kT}-1}\frac{R_{\ast}^{2}}{r(x)^{2}}\delta\left(\theta-0\right)$ (7)

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

$\tau_{\gamma\gamma}(\epsilon_{\gamma},\chi)=\Phi(\epsilon_{\gamma})\int_{0}^{2\pi}{\rm
d}\phi\int_{0}^{l}{\rm
d}x\,\frac{1-\mu(x,\chi)}{r(x,\chi)^{2}}=2\pi\Phi(\epsilon_{\gamma})f(\chi),$
(8)

where

$\Phi(\epsilon_{\gamma})=\frac{64\pi}{3h^{3}c^{3}}\sigma_{\rm
T}\left(\frac{m_{\rm
e}c^{2}}{\epsilon_{\gamma}}\right)^{3}\frac{R_{\ast}^{2}}{e^{2m_{\rm
e}c^{2}/\epsilon_{\gamma}kT}-1}$ (9)

is only a function of the gamma-ray energy and is constant with the orbital
phase. The  function

$f(\chi)=\int_{0}^{l}{\rm d}x\,\frac{1-\mu(x,\chi)}{r(x,\chi)^{2}}$ (10)

reflects the geometry of the system and determines the variability of the
optical depth with the orbital phase. The terms $\mu(x,\chi)$ and $r(x,\chi)$
can be expressed as functions of the distance between the star and gamma-ray
emitting region $d$ and the initial interaction angle,
$\mu_{0}\equiv\mu(x=0)$, as

$\mu(x,\chi)=\sqrt{\left(1-\frac{d(\chi)^{2}(1-\mu_{0}(\chi)^{2})}{r(x,\chi)^{2}}\right)},$
(11)

and

$r(x,\chi)^{2}=d(\chi)^{2}+x^{2}+2xd(\chi)\mu_{0}(\chi).$ (12)

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(\chi)$ in Eq. 10 and performing the integration we find that
$f(\chi)$ can be expressed analytically as (see Appendix A for details)

$f(\chi)=\frac{1}{d(\chi)}\left[\frac{1}{\sqrt{1-\mu_{0}(\chi)^{2}}}\left(\frac{\pi}{2}-\arctan{\frac{\mu_{0}(\chi)}{\sqrt{1-\mu_{0}(\chi)^{2}}}}\right)-1\right].$
(13)

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(\chi)$ and
$\mu_{0}(\chi)$ in terms of the orbital parameters of the system:

$\mu_{0}(\chi)=\sin{i}\sin{(\chi+\omega)},$ (14)
$d(\chi)=\frac{a(1-e^{2})}{1-e\cos{\chi}}.$ (15)

Following these simple considerations it can be seen that the maximum of the
function $f(\chi,i,\omega,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.

###  2.3 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
$\gamma\gamma$ optical depth $\tau_{\gamma\gamma}$ as a function of the
orbital phase $\chi$ with color coding indicating the range of varied
parameters: left panels - varied inclination with the longitude of periastron
fixed at $135^{\circ}$ and eccentricity fixed at $0.5$; middle panels - varied
longitude of periastron with the inclination fixed at $45^{\circ}$ and the
eccentricity fixed at $0.5$; and right panels - varied eccentricity with the
inclination fixed at $45^{\circ}$ and the longitude of periastron fixed at
$135^{\circ}$. 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
$\omega\sim 270^{\circ}$, 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 $\omega\sim 270^{\circ}$ (Fig. 3).

Figure 2: Optical depth $\tau_{\gamma\gamma}$ 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: Same as the right bottom plot in Fig. 2 but for
$\omega=270^{\circ}$.  Table 1: Physical parameters of examined systems.

| Test Binary | PSR B1259-63✫ | LS 5039✦ | PSR J2032$+$4127✝
---|---|---|---|---
$i$ [∘] | 45 | 22.2 | 60 | $>27.7$
$\omega$ [∘] ✠ | 135 | 138.7 | 45.8 (56.0) | $21-52$
$e$ | 0.5 | 0.87 | 0.35 (0.337) | $0.936-0.989$
$P$ [days] | 100 | $1,236.72$ | 3.906 | $16,000-17,670$
$T_{\ast}$ [K] | 30,000 | 33,000 | 39,000 | 20,000
$M_{\ast}$ [M⊙] | 30 | 31 | 22.9 | 15
$R_{\ast}$ [R⊙] | 10 | 9.2 | 9.3 | 10

✠ Longitude of periasrton of the pulsar orbit
✫ Orbital and stellar parameters are adopted from Negueruela et al. (2011) and
references therein
✦ Orbital and stellar parameters are adopted from Casares et al. (2005);
alternative orbital parameters (in brackets) are from Aragona et al. (2009)
✝ Ho et al. (2017) for orbital parameters and Lyne et al. (2015) for stellar
parameters

Figure 4: 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
B$1259-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.

Figure 4 shows the application of the optical depth calculation to two well
studied binary systems, PSR B$1259-63$ (top panel) and LS 5039 (bottom panel),
with substantially different orbital parameters (Table 1). While PSR
B$1259-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.

###  2.4 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 $\omega$ we calculate the orbital
phase with the highest optical depth, i.e. maximum absorption,
$\phi_{\mathrm{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^{\circ}$ on the estimate of $\phi_{\mathrm{max}}$.

##  3 Application to PSR J2032+4127/MT91 213

PSR J$2032+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 J$2032+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 $\sim
15M_{\odot}$ Be star MT$91213$ (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
($e>0.93$) and long  orbital period of about $45-50$ years (Ho et al., 2017).
The binary benefited from an extensive coverage across the electromagnetic
regime during its periastron passage (Abeysekara et al., 2018; Coe et al.,
2019; Ng et al., 2019; Coe et al., 2019), which undoubtedly provided a better
understanding of the system and allowed some constraints to be set on the
physical parameters. However, PSR J$2032+4127$/MT91 213 still suffers from a
poorly constrained orbit. The only available orbital solution provided by Ho
et al. (2017) contains a lot of uncertainties refining the orbital period in
the range $16000-17670$ days, eccentricity in the range $0.936-0.989$,
longitude of periastron in the range of $21^{\circ}-52^{\circ}$ and the
inclination angle above $27.7^{\circ}$ (taking into account the uncertainty on
the mass function and on the mass of the Be star).

VHE observations of the PSR J$2032+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 $\sim 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 X-ray 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 gamma-gamma absorption scenario for the
attenuation of the TeV flux.

To probe the orbital parameters of PSR J$2032+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 $\omega$. 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
$\phi_{\mathrm{max}}$ for higher eccentricity. The estimate of
$\phi_{\mathrm{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.

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).

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 J$2032+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^{\circ}-290^{\circ}$, 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 gamma-gamma absorption, i.e. those that fall
into the $186^{\circ}-290^{\circ}$ 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 $\omega$ 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-\omega$ 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.

Figure 6: Allowed combinations of values of $i$ and $\omega$ as derived from
the position of the minimum in the TeV light curve of PSR J$2032+4127$/MT91
213 that takes place at $7\pm 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 $\omega$ set by optical observations Ho
et al. (2017)

.

##  4 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 J$2032+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.

Iurii Sushch acknowledges support by the National Research Foundation of South
Africa (Grant Number 132276). Brian van Soelen acknowledges support by the
National Research Foundation of South Africa (Grant Number 119430). Some of
the numerical calculations presented were performed using the University of
the Free State High Performance Computing Unit.

##  Appendix A Derivation of function $f(\phi)$

Substituting Eqs. 11 and 12 $f(\phi)$ in Eq. 8 and changing the upper limit of
integration to infinity we get

$f=\int_{0}^{\infty}\frac{1-\frac{x+d\mu_{0}}{\sqrt{d^{2}+x^{2}+2xd\mu_{0}}}}{d^{2}+x^{2}+2xd\mu_{0}}{\rm
d}x.$ (A1)

We further substitute $y=x+d\mu_{0}$ and split it into two integrals

$\displaystyle f$ $\displaystyle=\int_{d\mu_{0}}^{\infty}\frac{{\rm
d}y}{d^{2}(1-\mu_{0}^{2})+y^{2}}-\int_{d\mu_{0}}^{\infty}\frac{y{\rm
d}y}{(d^{2}(1-\mu_{0}^{2})+y^{2})^{3/2}}$ (A2) $\displaystyle=I_{1}-I_{2}$

For the first integral $I_{1}$,

$\displaystyle I_{1}$ $\displaystyle=\int_{d\mu_{0}}^{\infty}\frac{{\rm
d}y}{d^{2}(1-\mu_{0}^{2})+y^{2}}$
$\displaystyle=\left.\frac{1}{d\sqrt{1-\mu_{0}^{2}}}\arctan{\frac{y}{d\sqrt{1-\mu_{0}^{2}}}}\right\rvert_{d\mu_{0}}^{\infty}$
(A3)
$\displaystyle=\frac{1}{d\sqrt{1-\mu_{0}^{2}}}\left(\frac{\pi}{2}-\arctan{\frac{\mu_{0}}{\sqrt{1-\mu_{0}^{2}}}}\right).$

For the second integral $I_{2}$, conducting another change of variables,
$t=y^{2}+d^{2}(1-\mu_{0}^{2})$, we have

$\displaystyle I_{2}$
$\displaystyle=\frac{1}{2}\int_{d^{2}}^{\infty}t^{-3/2}{\rm
d}t=\left.-t^{-1/2}\right\rvert_{d^{2}}^{\infty}=\frac{1}{d}$ (A4)

Combining A and A4 we get Eq. 13.

##  Appendix 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 $\eta=\cos{\theta}$ and limited to
$\eta\in\left[\cos{\theta_{\ast}};1\right]$, with $\theta_{\ast}$ being the
angular radius of the star as viewed from the gamma-ray location, i.e.
$\cos{\theta_{\ast}}=\frac{d}{\sqrt{d^{2}+R_{\ast}^{2}}}$, where $d$ is the
distance between the gamma-ray location and the center of the star and
$R_{\ast}$ 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 $\phi$ from 0 to $\pi$ using the symmetry. The solid angle integral in
Eq. 2 can be then rewritten as

${\int_{4\pi}{\rm
d}\Omega=\int_{0}^{2\pi}d\phi\int_{0}^{\theta_{\ast}}\sin{\theta}d\theta=2\int_{0}^{\pi}d\phi\int_{\eta_{\ast}}^{1}d\eta}$
(B1)

where $\eta_{\ast}=\cos{\theta_{\ast}}$. 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.

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
B$1259-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.

##  Appendix C Comparison of analytic and numeric approach

To calculate the orbital phase with the maximum absorption,
$\phi_{\mathrm{max}}$, we split the orbit into 200 bins setting a systematic
error of $1.8^{\circ}$. We estimate $\phi_{\mathrm{max}}$ for different
combinations of $i$ and $\omega$ 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 $\phi_{\mathrm{max}}$
as a function of $i$ (y-axis) and $\omega$ (x-axis) calculated numerically
(top panels) and analytically (upper middle panels). The lower middle panels
show the difference between two estimates,
$\phi_{\mathrm{max}}^{\mathrm{num}}-\phi_{\mathrm{max}}^{\mathrm{an}}$, and
the bottom panels map the the maximum optical depth as obtained from numerical
simulations, i.e. the optical depth at $\phi_{\mathrm{max}}^{\mathrm{num}}$.
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 $\phi_{\mathrm{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) $\omega$ is
irrelevant for the estimate of $\phi_{\mathrm{max}}$ as it stays at $\approx
180^{\circ}$ (periastron) determined by the  to the smallest binary separation
at this phase. For $i\gtrsim 10^{\circ}$, $\phi_{\mathrm{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 $\omega$ is close to $270^{\circ}$ (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 $\omega\sim 270^{\circ}$
features two maxima that further complicates things.

Apart for a small domain around $\omega\simeq 270^{\circ}$ 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 $\omega\simeq 270^{\circ}$ the analytic
solution becomes significantly different (although still at most by $\sim
15^{\circ}$) 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
$\tau_{\gamma\gamma}(\phi_{\mathrm{max}}^{\mathrm{num}})$. 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 $\phi_{\mathrm{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.

Figure 8: The upper two rows show the maps of the orbital phase of maximum
absorption, $phi_{\mathrm{max}}$, color-coded from $0^{\circ}$ to
$360^{\circ}$, 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,
$\phi_{\mathrm{max}}^{\mathrm{num}}-\phi_{\mathrm{max}}^{\mathrm{an}}$, and
the bottom maps show numerically calculated maximum optical depth,
$\tau_{\gamma\gamma}(\phi_{\mathrm{max}}^{\mathrm{num}})$. Figure 9: Same as
Fig. 8 but for different values of the orbital period and $e=0.5$

## References

  * Abdo et al. (2009) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, Science, 325, 840, doi: 10.1126/science.1175558
  * Abeysekara et al. (2018) Abeysekara, A. U., Benbow, W., Bird, R., et al. 2018, ApJ, 867, L19, doi: 10.3847/2041-8213/aae70e
  * Aharonian et al. (2005) Aharonian, F., Akhperjanian, A. G., Aye, K.-M., et al. 2005, A&A, 442, 1, doi: 10.1051/0004-6361:20052983
  * Aharonian et al. (2009) Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009, A&A, 507, 389, doi: 10.1051/0004-6361/200912339
  * Aragona et al. (2009) Aragona, C., McSwain, M. V., Grundstrom, E. D., et al. 2009, ApJ, 698, 514, doi: 10.1088/0004-637X/698/1/514
  * Boettcher et al. (2012) Boettcher, M., Harris, D. E., & Krawczynski, H. 2012, Relativistic Jets from Active Galactic Nuclei (Wiley)
  * Bosch-Ramon et al. (2015) Bosch-Ramon, V., Barkov, M. V., & Perucho, M. 2015, A&A, 577, A89, doi: 10.1051/0004-6361/201425228
  * Bosch-Ramon & Khangulyan (2009) Bosch-Ramon, V., & Khangulyan, D. 2009, International Journal of Modern Physics D, 18, 347, doi: 10.1142/S0218271809014601
  * Böttcher & Dermer (2005) Böttcher, M., & Dermer, C. D. 2005, ApJ, 634, L81, doi: 10.1086/498615
  * Camilo et al. (2009) Camilo, F., Ray, P. S., Ransom, S. M., et al. 2009, ApJ, 705, 1, doi: 10.1088/0004-637X/705/1/1
  * Casares et al. (2005) Casares, J., Ribó, M., Ribas, I., et al. 2005, MNRAS, 364, 899, doi: 10.1111/j.1365-2966.2005.09617.x
  * Chernyakova & Malyshev (2020) Chernyakova, M., & Malyshev, D. 2020, in Multifrequency Behaviour of High Energy Cosmic Sources - XIII. 3-8 June 2019. Palermo, 45, doi: 10.22323/1.362.0045
  * Chernyakova et al. (2019) Chernyakova, M., Malyshev, D., Paizis, A., et al. 2019, A&A, 631, A177, doi: 10.1051/0004-6361/201936501
  * Coe et al. (2019) Coe, M. J., Okazaki, A. T., Steele, I. A., et al. 2019, MNRAS, 485, 1864, doi: 10.1093/mnras/stz515
  * Dermer & Böttcher (2006) Dermer, C. D., & Böttcher, M. 2006, ApJ, 643, 1081, doi: 10.1086/502966
  * Dubus (2006a) Dubus, G. 2006a, A&A, 456, 801, doi: 10.1051/0004-6361:20054779
  * Dubus (2006b) —. 2006b, A&A, 451, 9, doi: 10.1051/0004-6361:20054233
  * Dubus (2013) —. 2013, A&A Rev., 21, 64, doi: 10.1007/s00159-013-0064-5
  * Dubus (2015) —. 2015, Comptes Rendus Physique, 16, 661, doi: 10.1016/j.crhy.2015.08.014
  * Dubus et al. (2015) Dubus, G., Lamberts, A., & Fromang, S. 2015, A&A, 581, A27, doi: 10.1051/0004-6361/201425394
  * Gould & Schréder (1967) Gould, R. J., & Schréder, G. P. 1967, Physical Review, 155, 1404, doi: 10.1103/PhysRev.155.1404
  * H. E. S. S. Collaboration et al. (2020) H. E. S. S. Collaboration, Abdalla, H., Adam, R., et al. 2020, A&A, 633, A102, doi: 10.1051/0004-6361/201936621
  * H.E.S.S. Collaboration et al. (2013) H.E.S.S. Collaboration, Abramowski, A., Acero, F., et al. 2013, A&A, 551, A94, doi: 10.1051/0004-6361/201220612
  * Ho et al. (2017) Ho, W. C. G., Ng, C.-Y., Lyne, A. G., et al. 2017, MNRAS, 464, 1211, doi: 10.1093/mnras/stw2420
  * Huber et al. (2021a) Huber, D., Kissmann, R., Reimer, A., & Reimer, O. 2021a, A&A, 646, A91, doi: 10.1051/0004-6361/202039277
  * Huber et al. (2021b) Huber, D., Kissmann, R., & Reimer, O. 2021b, A&A, 649, A71, doi: 10.1051/0004-6361/202039278
  * Jauch & Rohrlich (1976) Jauch, J. M., & Rohrlich, F. 1976, The theory of photons and electrons. The relativistic quantum field theory of charged particles with spin one-half
  * Khangulyan et al. (2008) Khangulyan, D., Aharonian, F., & Bosch-Ramon, V. 2008, MNRAS, 383, 467, doi: 10.1111/j.1365-2966.2007.12572.x
  * Lyne et al. (2015) Lyne, A. G., Stappers, B. W., Keith, M. J., et al. 2015, MNRAS, 451, 581, doi: 10.1093/mnras/stv236
  * Negueruela et al. (2011) Negueruela, I., Ribó, M., Herrero, A., et al. 2011, ApJ, 732, L11, doi: 10.1088/2041-8205/732/1/L11
  * Ng et al. (2019) Ng, C. Y., Ho, W. C. G., Gotthelf, E. V., et al. 2019, ApJ, 880, 147, doi: 10.3847/1538-4357/ab2adb
  * Olivera Nieto (2023) Olivera Nieto, L. 2023, PhD thesis, MPIK, Heidelberg, doi: 10.11588/heidok.00032936
  * Romoli et al. (2015) Romoli, C., Bordas, P., Mariaud, C., et al. 2015, ArXiv e-prints.  https://arxiv.org/abs/1509.03090
  * Sushch & van Soelen (2017) Sushch, I., & van Soelen, B. 2017, ApJ, 837, 175, doi: 10.3847/1538-4357/aa62ff
  * Weng et al. (2022) Weng, S.-S., Qian, L., Wang, B.-J., et al. 2022, Nature Astronomy, 6, 698, doi: 10.1038/s41550-022-01630-1
  * Zdziarski & Lightman (1985) Zdziarski, A. A., & Lightman, A. P. 1985, ApJ, 294, L79, doi: 10.1086/184513


# On Imperfect Recall in Multi-Agent Influence Diagrams

James Fox University of Oxford<EMAIL_ADDRESS>Imperial College
LondonUniversity of OxfordUniversity of OxfordUniversity of OxfordUniversity
of Oxford Matt MacDermott Imperial College London
<EMAIL_ADDRESS>University of OxfordUniversity of
OxfordUniversity of OxfordUniversity of Oxford Lewis Hammond University of
Oxford<EMAIL_ADDRESS>University of OxfordUniversity of
OxfordUniversity of Oxford Paul Harrenstein University of Oxford
<EMAIL_ADDRESS>University of OxfordUniversity of Oxford
Alessandro Abate University of Oxford<EMAIL_ADDRESS>University of Oxford
Michael Wooldridge University of Oxford<EMAIL_ADDRESS>
###### Abstract

Multi-agent influence diagrams (MAIDs) are a popular game-theoretic model
based on Bayesian networks. In some settings, MAIDs offer significant
advantages over extensive-form game representations. Previous work on MAIDs
has assumed that agents employ behavioural policies, which set independent
conditional probability distributions over actions for each of their
decisions. In settings with imperfect recall, however, a Nash equilibrium in
behavioural policies may not exist. We overcome this by showing how to solve
MAIDs with forgetful and absent-minded agents using mixed policies and two
types of correlated equilibrium. We also analyse the computational complexity
of key decision problems in MAIDs, and explore tractable cases. Finally, we
describe applications of MAIDs to Markov games and team situations, where
imperfect recall is often unavoidable.

##  1 Introduction

Multi-agent influence diagrams (MAIDs) are a graphical representation for
dynamic non-cooperative games, which can be more compact and expressive than
extensive-form games (EFGs) [koller2003multi]. Like Bayesian networks (BNs),
MAIDs use a directed acyclic graph (DAG) to represent conditional
probabilistic dependencies between random variables, but they also specify
decision and utility variables for each agent. Each agent selects a
behavioural policy – independent conditional probability distributions (CPDs)
over actions for each of their decision variables – to maximise their expected
utility. A MAID’s mechanised graph extends this DAG by explicitly representing
each variable’s distribution and showing which other variables’ distributions
matter to an agent optimising a particular decision rule [causalgames,
koller2003multi, dawid2002influence].

MAIDs, and their causal variants [causalgames], have been used in the design
of safe and fair AI systems [everitt2021agent, ashurst2022fair,
everitt2021reward, farquhar2022path, carroll2023characterizing], to explore
reasoning patterns and deception [pfeffer2007reasoning, ward2022agent], and to
identify agents from data [kenton2022discovering]. However, to date, agents in
MAIDs are usually assumed to have perfect (or, at least, ‘sufficient’) recall
[koller2003multi]. This assumption is often unreasonable. For example, MAIDs
must allow imperfect recall to handle bounded rationality, teams with
imperfect communication [detwarasiti2005influence], or memoryless policies in
Markov games. However, forgetfulness (of previous observations) or absent-
mindedness (about whether previous decisions have even been made) can prevent
the existence of a Nash Equilibrium (NE) in behavioural policies. To overcome
this, one can consider other solution concepts, such as mixed or correlated
equilibria.

In this work, we focus on imperfect recall in MAIDs. Imperfect recall has
already been extensively studied in EFGs [piccione1997interpretation,
Kuhn1953, waugh2009practical], but a MAID’s mechanised graph makes graphically
explicit the semantic difference between behavioural and mixed policies
(hidden in EFGs) and readily identifies forgetful or absent-minded agents (or
teams). Our insights inspire two definitions of _correlated equilibrium_ in
MAIDs. The first follows from the normal-form game definition
[aumann1974subjectivity]. The second, based on von Stengel and Forges’
extensive-form correlated equilibrium [von2008extensive], is more natural for
dynamic settings, can yield greater social welfare, and is easier to compute.
Again, mechanised graphs clearly depict the assumptions made in both. Next, we
examine MAIDs from a computational complexity perspective by studying the
decision problems of finding a best response, checking whether a policy
profile is an NE, and checking whether each type of NE exists. These provide
an insight into what makes particular instances hard, when computations can be
made tractable, and rigorously identify which problems are suitable for
analysis as MAIDs. Our results also apply to refinements of MAIDs, such as
causal games [causalgames]. We assume familiarity with EFGs
[maschler2020game], BNs [koller2009probabilistic], and the complexity classes
$\mathrm{P}$, $\mathrm{NP}$, and $\mathrm{PP}$ [papadimitriou1994]. Proof
sketches are provided, but details are deferred to the appendices.

#### Related Work.

There is a rich literature on influence diagrams [kjaerulff2008bayesian] and
imperfect recall has been studied in single-agent influence diagrams
[maua2012solving, maua2016fast, lauritzen2001representing, de2008strategy,
van2022complete] as well as in EFGs [aumann1997absent, kaneko1995behavior,
Kuhn1953, piccione1997interpretation, waugh2009practical]. However, to our
knowledge, we are the first to focus on imperfect recall in influence diagrams
with multiple agents.

A full policy profile in a MAID induces a BN, so many of our results inherit
from that setting, where the decision problem variant of marginal inference
is, in general, $\mathrm{PP}$-complete [littman2001stochastic]. However, we
care about the cases we encounter in practice, not just the worst case.
Marginal inference in a BN can be performed in time exponential in the
treewidth of the underlying graph [koller2009probabilistic], which entails a
poly-time algorithm when the treewidth is small. Similarly, we will see that
tractable results for computations in MAIDs can be found when problems are
restricted to certain settings. We also sometimes reduce from partial order
games [zahoransky2021partial], which can be interpreted as MAIDs without
chance nodes, with deterministic decision rules, and where each agent has a
single utility node as a child of all the decision nodes.

##  2 The Model

We use capital letters $V$ for random variables, lowercase letters $v$ for
their instantiations, and bold letters $\bm{V}$ and $\bm{v}$, respectively,
for sets of variables and their instantiations. We let $\textit{dom}(V)$
denote the (finite, non-singleton) domain of $V$ (for ease, we take this to be
binary unless stated otherwise) and
$\textit{dom}(\bm{V})\coloneqq\bigtimes_{V\in\bm{V}}\textit{dom}(V)$. Parents
and children of $V$ in a graph are denoted by $\textbf{Pa}_{V}$ and
$\textbf{Ch}_{V}$, respectively (with $\textbf{pa}_{V}$ and $\textbf{ch}_{V}$
their instantiations) and $\Delta(X)$ denotes the set of all probability
distributions over a set $X$.

######  Example 1.

An autonomous taxi decides whether to offer Alice a discount ($T$) depending
on whether its journey count exceeds a quota ($Q$). Alice decides whether to
accept a journey ($A$) depending on the price. The taxi wants to maximise
profit, but if its journey count is less than the quota and Alice rejects it,
the taxi pays a penalty (the municipality uses this mechanism to prevent a
proliferation of unnecessary taxis). Alice’s utility is a function of her
decision and the price offered by the taxi.

Figure 1(a) shows a MAID for this example. Chance variables (moves by nature),
decision variables, and utility variables are represented by white circles,
squares, and diamonds, respectively. Full edges leading into chance and
utility nodes represent probabilistic dependence, as in a BN. Dotted edges
leading into decision nodes identify information available to the agent when a
decision $D$ is made, so $\textbf{pa}_{D}$, the values of $\textbf{Pa}_{D}$,
represents the decision context for $D$. In EFGs, imperfect information is
represented using explicitly labelled information sets. In MAIDs, we can infer
that Alice is unaware of the value of $Q$ when making her decision by the lack
of edge $Q\rightarrow A$. A parameterisation defines the CPDs for the chance
and utility variables, whereas CPDs of decision nodes are chosen by the agents
playing the game.

######  Definition 1 ([koller2003multi]).

A multi-agent influence diagram (MAID) is a structure
$\mathcal{M}=(\mathcal{G},\bm{\theta})$. $\mathcal{G}=(N,\bm{V},E)$ specifies
a set of agents $N=\\{1,\dots,n\\}$ and a DAG $(\bm{V},E)$, where $\bm{V}$ is
partitioned into chance variables $\bm{X}$, decision variables
$\bm{D}=\bigcup_{i\in N}\bm{D}^{i}$, and utility variables
$\bm{U}=\bigcup_{i\in N}\bm{U}^{i}$. The parameters
$\bm{\theta}=\\{\theta_{V}\\}_{V\in\bm{V}\setminus\bm{D}}$ define the CPDs
$\Pr(V\mid\textbf{Pa}_{V})$ for each non-decision variable such that for any
setting of the decision variables’ CPDs, the resulting joint distribution over
$\bm{V}$ is Markov compatible with the DAG, i.e.,
$\Pr(\bm{v})=\prod_{V\in\bm{V}}\Pr(v\mid\textbf{pa}_{V})$.

$Q$$T$$A$$U^{T}$$U^{A}$Bernoulli(0.5)$U^{T}=T\cdot
A-J(1-A)$$U^{A}=A\cdot(3-T)$

(a)

$Q$$T$$A$$\Theta_{Q}$$\Pi_{T}$$\Pi_{A}$$U^{T}$$U^{A}$$\Theta_{U^{T}}$$\Theta_{U^{A}}$

(b)

Figure 1: A MAID (a) and its mechanised graph (b) for Example 1, which is a
perfect recall and imperfect, but sufficient, information game.

Given a MAID, a decision rule $\pi_{D}$ for $D\in\bm{D}$ is a CPD
$\pi_{D}(D\mid\textbf{Pa}_{D})$. A partial (behavioural) policy profile
$\pi_{\bm{D}^{\prime}}$ is a set of decision rules for each
$D\in\bm{D}^{\prime}\subseteq\bm{D}$, whereas $\pi_{-\bm{D}^{\prime}}$ is the
set of decision rules for each $D\in\bm{D}\setminus\bm{D}^{\prime}$. A
(behavioural) policy ${\bm{\pi}}^{i}$ refers to ${\bm{\pi}}_{\bm{D^{i}}}$, and
a (full) policy profile ${\bm{\pi}}=({\bm{\pi}}^{1},\ldots,{\bm{\pi}}^{n})$ is
a tuple of policies, where
${\bm{\pi}}^{-i}\coloneqq({\bm{\pi}}^{1},\dots,{\bm{\pi}}^{i-1},{\bm{\pi}}^{i+1},\dots,{\bm{\pi}}^{n})$.
A decision rule is pure if $\pi_{D}(d\mid\textbf{pa}_{D})\in\\{0,1\\}$, which
holds for a policy (profile) if it holds for all decision rules in the policy
(profile). For clarity, we use an overhead dot to mark this determinism, e.g.,
$\dot{\pi}_{D},\dot{{\bm{\pi}}}^{i}$, or $\dot{{\bm{\pi}}}$.

By combining ${\bm{\pi}}$ with the partial distribution $\Pr$ over the chance
and utility variables, we obtain a joint distribution:

$\Pr^{\bm{\pi}}(\bm{x},\bm{d},\bm{u})\coloneqq\prod_{V\in\bm{V}\setminus\bm{D}}\Pr(v\mid\textbf{pa}_{V})\cdot\prod_{D\in\bm{D}}\pi_{D}(d\mid\textbf{pa}_{D})$

A full policy profile ${\bm{\pi}}$ therefore induces a BN with DAG given by
the MAID’s graph. Agent $i$’s expected utility $EU^{i}({\bm{\pi}})$ for a
given policy profile ${\bm{\pi}}$ is defined as the expected sum of their
utility variables:

$\textstyle
EU^{i}({\bm{\pi}})\coloneqq\sum_{U\in\bm{U}^{i}}\sum_{u\in\textit{dom}(U)}\Pr^{{\bm{\pi}}}(U=u)\cdot
u$

Utility variables have deterministic CPDs, so can be interpreted as functions
$U:\textit{dom}(\textbf{Pa}_{U})\rightarrow~{}\mathbb{R}$ to show their
functional dependence on their parents (e.g., Figure 1(a)). An NE is defined
in the usual way.

######  Definition 2 ([koller2003multi]).

A (behavioural) policy profile ${\bm{\pi}}$ is a Nash equilibrium (NE) (in
behavioural policies) if for every agent $i\in N$ and every alternative
(behavioural) policy ${\bm{\varpi}}^{i}$:
$EU^{i}({\bm{\pi}^{-i},{\bm{\pi}}^{i}})\geq
EU^{i}({\bm{\pi}^{-i},{\bm{\varpi}}^{i}})$

Collectively, the decision rules of decision variables and the CPDs of chance
or utility nodes are known as mechanisms. A mechanism $\mathsf{M}_{V}$ for $V$
is strategically relevant to a decision rule for $D$ if the choice of the CPD
at $\mathsf{M}_{V}$ can affect the optimal choice of this decision rule.
Koller and Milch [koller2003multi] define an associated sound and complete
graphical criterion for strategic relevance, $\bm{s}$-reachability, based on
d-separation which can be checked in $\mathcal{O}(|\bm{V}|+|E|)$ time
[shachter2013bayes] (see Appendix LABEL:app:stratrel for formal definitions).

A MAID’s regular graph $\mathcal{G}$ captures the probabilistic dependencies
between object-level variables in the game’s environment, but its mechanised
graph $\mathsf{m}{\mathcal{G}}$ is an enhanced representation which adds an
explicit representation of the strategically relevant dependencies between
agents’ decision rules and the game’s parameterisation (see [causalgames] for
details). Each object-level variable $V\in\bm{V}$ has a mechanism parent
$\mathsf{M}_{V}$ representing the distribution governing $V$: each decision
$D$ has a new _decision rule_ parent $\Pi_{D}=\mathsf{M}_{D}$ and each non-
decision $V$ has a new _parameter_ parent $\Theta_{V}=\mathsf{M}_{V}$, whose
values parameterise the CPDs.

Agents select a decision rule $\pi_{D}$ (i.e., the value of a decision rule
variable $\Pi_{D}$) based on both the parameterisation of the game (i.e., the
values of the parameter variables) and the selection of the other decision
rules ${\bm{\pi}}_{-D}$ – these dependencies are captured by the edges from
other mechanisms into decision rule nodes. $s$-reachability determines which
of these edges are necessary, so $\mathsf{M}_{V}\rightarrow\Pi_{D}$ exists if
and only if $\Pi_{D}$ strategically relies on $\mathsf{M}_{V}$. The mechanised
graph for Example 1 (in Figure 1(b)) shows that $\Pi_{T}$ strategically relies
on $\Theta_{U^{T}}$ and $\Pi_{A}$, whereas $\Pi_{A}$ only strategically relies
on $\Theta_{U^{A}}$. In contrast to a MAID’s regular graph $\mathcal{G}$,
which is a DAG, there may exist cycles between mechanisms (e.g., Figure
LABEL:fig:noNE:a).

For convenience, we denote the set of agent $i$’s behavioural policies as
${\bm{P}}^{i}\coloneqq\textit{dom}(\bm{\Pi}^{i})$, with sets of pure policies
denoted as $\dot{{\bm{P}}}^{i}$ and (pure) policy profiles denoted by
${\bm{P}}$ ($\dot{{\bm{P}}}$).

###  2.1 Concise Representations

A concise representation of MAIDs is needed for three reasons. First, real
numbers may obscure the true complexity of the problems
[bodlaender2002complexity], so we assume that all probability parameters are
given by a fraction of two integers, both expressed in finite binary notation.
This is realistic since the probabilities are normally either assessed by
domain experts or estimated by a learning algorithm and means that all CPDs
can be read in poly-time. Second, even with binary variables, a joint
distribution across $\bm{V}$ requires $2^{|\bm{V}|}-1$ parameters. A MAID or
BN’s graphical Markov factorisation reduces this to
$\sum_{V\in\bm{V}}2^{|\textbf{Pa}_{V}|}$, but this can still be exponential in
$|\bm{V}|$. Therefore, it is standard [shimony1994finding, roth1996hardness,
kwisthout2009computational, koller2009probabilistic] to assume that the
maximum in-degree in the graph is much less than $|\bm{V}|$ (or constant), so
that the size of the CPDs are polynomial in $|\bm{V}|$. This means that the
total representation of our MAID (including all CPDs) is polynomial in our
chosen complexity parameter $|\bm{V}|$. Finally, as in BNs, our complexity
results are strongly affected by the DAG’s treewidth. The treewidth of a DAG
measures its resemblance to a tree and is given by the number of vertices in
the largest clique of the corresponding triangulated moral graph minus one
[bodlaender1993linear].

##  3 Imperfect Recall in MAIDs

Agents may possess different degrees of information about the state of a game.
A game has perfect recall if each agent remembers all their past decisions and
observations, and it has perfect information if each agent is aware of every
agent’s past decisions and observations.

######  Definition 3 ([koller2003multi]).

Agent $i$ in a MAID $\mathcal{M}$ is said to have perfect recall if there
exists a total ordering $D_{1}\prec\cdots\prec D_{m}$ over $\bm{D}^{i}$ such
that $(\textbf{Pa}_{D_{j}}\cup D_{j})\subseteq\textbf{Pa}_{D_{k}}$ for any
$1\leq j<k\leq m$. $\mathcal{M}$ is a perfect recall game if all agents in
$\mathcal{M}$ have perfect recall. $\mathcal{M}$ is a perfect information game
if there exists such an ordering over $\bm{D}$.

A MAID with perfect information (recall) can be transformed into an EFG with
perfect information (recall), and vice versa [Hammond2021]. Hence, these
information conditions also guarantee the existence of an NE in pure
(behavioural) policies in the MAID ([Kuhn1953] gives the equivalent results in
EFGs). However, the mechanised representation of a MAID enables weaker
criteria to be defined – sufficient information and sufficient recall. Later,
in Proposition LABEL:prop:suffexistence, we will see that these criteria
preserve the NE existence results of perfect information and perfect recall
games, respectively.

######  Definition 4.

Agent $i$ in a MAID $\mathcal{M}$ has sufficient recall [milch2008ignorable]
if the subgraph of the mechanised graph $\mathsf{m}{\mathcal{G}}$ restricted
to just agent $i$’s decision rule nodes $\bm{\Pi}_{\bm{D}^{i}}$ is acyclic.
$\mathcal{M}$ is a sufficient recall game if all agents in $\mathcal{M}$ have
sufficient recall. $\mathcal{M}$ is a sufficient information game if the
subgraph of $\mathsf{m}{\mathcal{G}}$ restricted to contain only and all
decision rule nodes $\bm{\Pi}_{\bm{D}}$ is acyclic.111Note that since previous
work on influence diagrams has not modelled absent-mindedness (see our
Definition LABEL:def:forgetabsent in Section 3.1), this definition implicitly
assumes each mechanism variable has a single child.

###  3.1 Forgetfulness and Absent-Mindedness