11institutetext:  Department of Astronomy, School of Physics, Huazhong
University of Science and Technology, Wuhan 430074, China.
<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Modelling the correlated keV/TeV light curves of Be/$\gamma$-ray binaries

A. M. Chen 11 J. Takata 11

(Received xxx / Accepted xxx)

Be/$\gamma$-ray binaries comprise a confirmed or presumptive pulsar orbiting a
Be star and emit luminous $\gamma$-rays. Non-thermal emissions are thought to
arise from synchrotron radiation and inverse-Compton (IC) scattering in the
shock where the pulsar wind is terminated by the stellar outflow. We study
wind interactions and shock radiations from such systems and show that the
bimodal structures observed in keV/TeV light curves are caused by enhanced
synchrotron radiation and IC scattering during disc passages. We use a simple
radiation model to reproduce orbital modulations of keV X-ray and TeV
$\gamma$-ray flux and compare with two confirmed pulsar/Be star binaries (i.e.
PSR B1259-63/LS 2883 and PSR J2032+4127/MT91 213), and two candidates (i.e.
HESS J0632+057 and LS I +61∘303). We find that the keV/TeV light curves of
former two binaries can be well explained by the inclined disc model, while
modelling the modulated emissions of latter two sources remains challenging
with current orbital solutions. Therefore, we propose alternative orbital
geometries for HESS J0632+057 and LS I +61∘303\. We estimate the positions and
inclination angles of Be discs by fitting correlated keV/TeV light curves. Our
results could be beneficial for future measurements of orbital parameters and
searches for radio pulsations from presumed pulsars.

###### Key Words.:

binaries: close – X-rays: binaries – Gamma rays: stars – pulsars: individual

##  1 Introduction

$\gamma$-ray binaries are a rare subclass of high-mass binary systems
harbouring a compact object in orbit with a type O/B star and radiate luminous
$\gamma$-rays with non-thermal energy spectra peaking beyond 1 MeV (see Dubus
2013, 2015; Lamberts 2016; van Soelen 2017; Paredes & Bordas 2019; Chernyakova
& Malyshev 2020 for reviews). So far, less than ten such binaries have been
found, and only two of them with compact objects have been identified as
pulsars. The lack of accretion emission in other binaries indicates that the
unknown compact objects are likely to be non-accreting neutron stars (NSs),
although stellar-mass black holes cannot be completely ruled out (Dubus
2006a). Depending on the spectral types of massive companions and the presence
of a decretion disc, the detected $\gamma$-ray binaries can be divided into
the following two types:

  * •

O stars: such as LS 5039 (Motch et al. 1997; Aharonian et al. 2005a), 1FGL
J1018.6-5856 (Fermi LAT Collaboration 2012; H.E.S.S. Collaboration 2015a), LMC
P3 (Corbet et al. 2016; H.E.S.S. Collaboration 2018), HESS J1832-093 (H.E.S.S.
Collaboration 2015b; Eger et al. 2016), and 4FGL J1405.1-6119 (Corbet et al.
2019);

  * •

Be stars: such as PSR B1259-63/LS 2883111At first, the massive star, LS 2883,
was identified as a B2e-type star by Johnston et al. (1994), but it was
revised as type O9.5Ve by Negueruela et al. (2011). For simplicity, we still
attribute PSR B1259-63/LS 2883 as a Be/$\gamma$-ray binary. (Johnston et al.
1992; Aharonian et al. 2005b), PSR J2032+4127/MT91 213 (Lyne et al. 2015;
Abeysekara et al. 2018), HESS J0632+057/MWC 148 (Aharonian et al. 2007; Hinton
et al. 2009), and LS I +61∘303 (Gregory & Taylor 1978; Albert et al. 2006).

Both types show similar spectral characteristics with orbital modulated
radiations from radio to $\gamma$-rays, suggesting that they are likely
powered by spin-down of pulsars, as in PSR B1259-63/LS 2883. However, there
are significant differences between these two types. For example, those
$\gamma$-ray binaries hosting O stars usually have shorter orbital periods and
lower eccentricities than those with Be companions (hereafter Be/$\gamma$-ray
binaries). The former type usually exhibits a single-peak profile in its light
curves, which is typically attributed to Doppler-boosted shock emission around
inferior conjunction (INFC, Dubus et al. 2010; Takata et al. 2014; Molina &
Bosch-Ramon 2020). The keV/TeV light curves of Be/$\gamma$-ray binaries
usually display bimodal structures, which are likely caused by the compact
objects interacting with the discs of their massive companions (Tavani & Arons
1997; Chen et al. 2019).

Among detected Be/$\gamma$-ray binaries, PSR B1259-63/LS 2883 has been
thoroughly investigated both observationally and theoretically. The most
notable feature of this binary is the two-peak structures displayed in its keV
and TeV light curves. The double-hump behaviours associated with the
disappearance of radio pulsation around periastron are attributed to the
pulsar crossing the inclined Be disc (Chernyakova et al. 2006, 2014; Chen et
al. 2019, 2021a). The widely accepted scenario for PSR B1259-63/LS 2883
involves a pulsar wind shock terminated by the stellar outflow, where non-
thermal emissions are produced by the shock-accelerated particles (Kirk et al.
1999; Takata & Taam 2009; Kong et al. 2011; Takata et al. 2012). Chen et al.
(2019) proposed that the energy densities of magnetic field and photon field
would be enhanced during disc passages, increasing synchrotron and IC
luminosities, and resulting in the double peaks seen in the keV and TeV light
curves. Analogous behaviours are also seen in another pulsar/Be binary, PSR
J2032+4127/MT91 213 (Chernyakova et al. 2020b). The similar nature of these
two binaries suggests that a common physical mechanism (i.e. the inclined disc
model) is causing the orbital modulations of the multi-wavelength emissions.

The correlations and double-hump behaviours in the keV and TeV flux have been
found in another two Be/$\gamma$-ray binaries (i.e. HESS J0632+057 and LS I
+61∘303). However, the unknown natures of their compact objects and the lack
of well-measured orbital parameters make it challenging to understand the
origins and orbital modulations of their high-energy radiation. Several
emission models have been proposed, including: (1) the microquasar model,
where a pair of relativistic jets produce the non-thermal emission (e.g.
Bosch-Ramon & Paredes 2004; Romero et al. 2007; Zhang et al. 2009; Massi &
Torricelli-Ciamponi 2014); (2) the flip-flop scenario, where the compact
object is an NS that evolves from a rotational regime into a propeller state
along with the orbital motion (e.g. Torres et al. 2012; Papitto et al. 2012);
(3) the pulsar model, which involves a termination shock formed by collisions
between the pulsar wind and the stellar outflow as mentioned above (e.g.
Cerutti et al. 2008; Sierpowska & Torres 2009; Zdziarski et al. 2010; Zabalza
et al. 2011a; Bosch-Ramon et al. 2017; Malyshev et al. 2019).

Bimodal structures are seen in the keV/TeV light curves of all four discovered
Be/$\gamma$-ray binaries. In this paper, we attempt to study the origin of
their correlated keV/TeV emissions following the pulsar wind interaction model
and to investigate the formation of double-peak behaviours. The bimodal
structures of light curves could also provide indications as to the positions
and inclinations of Be discs in binaries, which are essential for measuring
the orbital parameters and searching for radio pulsations from the presumptive
pulsars. The paper is organised as follows: we describe the stellar outflow of
Be stars in Section 2.1; the termination shock and related radiation processes
are presented in Section 2.2 and 2.3, respectively. We then display our
calculated light curves with comparisons of observational data in Section 3.
Finally, a brief discussion and conclusions are presented in Section 4.

##  2 Model description

###  2.1 The stellar outflow

Figure 1: Illustration of the binary geometry with a pulsar orbiting around a
Be star.

Massive Be stars are characterised by Balmer emission lines and infrared (IR)
excess in their spectra (Porter & Rivinius 2003; Rivinius et al. 2013). The
origins of emission lines and IR excesses are attributed to a dense equatorial
disc expelled from the fast-rotating star. For Be/$\gamma$-ray binaries, the
stellar spin axis is usually misaligned with the orbital axis, and therefore
the equatorial disc is tilted relative to the orbital plane. In Fig. 1, we
illustrate the binary geometry of a pulsar in orbit with a Be companion. We
set the origin of Cartesian coordinates on the binary barycentre, the $x$-axis
towards periastron, and the $z$-axis along the orbital axis. Therefore, the
unit vectors in the directions of the pulsar, the observer, and the disc
normal can be written as follows (Chen et al. 2021a):

$\displaystyle\@vec{e}_{\rm{psr}}$ $\displaystyle=$
$\displaystyle\left(\cos\phi,\sin\phi,0\right),$ (1)
$\displaystyle\@vec{e}_{\rm{obs}}$ $\displaystyle=$ $\displaystyle\left(\sin
i_{\rm{o}}\cos\phi_{\rm{o}},\sin i_{\rm{o}}\sin\phi_{\rm{o}},\cos
i_{\rm{o}}\right),$ (2) $\displaystyle\@vec{e}_{\rm{disc}}$ $\displaystyle=$
$\displaystyle\left(\sin i_{\rm{d}}\cos\phi_{\rm{d}},\sin
i_{\rm{d}}\sin\phi_{\rm{d}},\cos i_{\rm{d}}\right),$ (3)

with $\phi$ being the pulsar’s true anomaly, and $i_{\rm o}$ ($i_{\rm d}$) and
$\phi_{\rm o}$ ($\phi_{\rm d}$) being the inclination angle and true anomaly
of the observer (the disc normal) projected on the orbital plane,
respectively.

Assuming that the stellar outflow is axisymmetric along the stellar spin axis,
we write the ram pressure of the outflow as (Ignace $\&$ Brimeyer 2006;
Petropoulou et al. 2018; Chen et al. 2019):

$\displaystyle p_{\rm w}(R,\theta)$ $\displaystyle=$ $\displaystyle
p_{0}R^{-2}(1+G\mid\cos\theta\mid^{m}),$ (4)

with $R$ being the radial distance and $\theta$ the latitude measured from the
equator. The value of $p_{0}$ is determined by the mass-loss rate $\dot{M}$
and velocity $v_{\rm w}$ of the polar wind as $p_{0}=\dot{M}v_{\rm w}$, and
the second term in the bracket depicts the additional pressure of the Be disc,
with $G$ being the pressure contrast. The half-opening angle of the Be disc is
related to the confinement parameter $m$ as (Ignace & Brimeyer 2006):

$\displaystyle\Delta\theta_{\rm d}$ $\displaystyle=$
$\displaystyle\arccos(2^{-1/m}),$ (5)

and the disc half-opening angle projected on the orbital plane is

$\displaystyle\Delta\phi_{\rm d}$ $\displaystyle=$
$\displaystyle\arcsin\left(\frac{\sin\Delta\theta_{\rm d}}{\sin{i_{\rm
d}}}\right).$ (6)

A well-confined disc requires a larger $m$, as depicted in the top panel of
Fig. 2. It is widely believed that the equatorial discs of Be stars are
geometrically thin, and so we simply adopt $m=100$ in calculations (which
corresponds to $\Delta\theta_{\rm d}=6.7^{\circ}$), and leave $i_{\rm d}$ as a
free model parameter. For a fixed value of $\Delta\theta_{\rm d}$, the disc
region projected on the orbital plane is larger with a smaller $i_{\rm d}$, as
presented in the bottom panel of Fig. 2. If $i_{\rm d}\leq\Delta\theta_{\rm
d}$, the disc is immersed on the plane, and then the pulsar will always be
interacting with the disc along the orbit. When $i_{\rm d}=90^{\circ}$, we
simply have $\Delta\phi_{\rm d}=\Delta\theta_{\rm d}$.

Figure 2: Top: Half-opening angle of the stellar disc $\Delta\theta_{\rm d}$
as a function of the confinement parameter $m$. Bottom: Half-opening angle of
the disc projected on the orbital plane $\Delta\phi_{\rm d}$ as a function of
the inclination angle $i_{\rm d}$, with the vertical and horizontal lines
corresponding to $\Delta\theta_{\rm d}=6.7^{\circ}$.

###  2.2 The termination shock

The energetic pulsar in a $\gamma$-ray binary drives a relativistic wind which
is terminated by stellar outflow. The position of the termination shock is
governed by the balance of two winds as:

$\displaystyle\frac{L_{\rm sd}}{4\pi r_{\rm s}^{2}c}$ $\displaystyle=$
$\displaystyle p_{0}R_{\rm s}^{-2}(1+G\mid\cos\theta\mid^{m}),$ (7)

where $L_{\rm sd}$ is the pulsar’s spin-down luminosity, and $c$ is the speed
of light. As the pulsar moves on the orbital plane, we have:

$\displaystyle\theta$ $\displaystyle=$
$\displaystyle\pi/2-\arccos\left(\@vec{e}_{\rm{disc}}\cdot\@vec{e}_{\rm{psr}}\right),$
(8)

and the disc midplane intersects on the orbital plane at true anomalies:

$\displaystyle\phi_{\rm d,\pm}$ $\displaystyle=$ $\displaystyle\phi_{\rm
d}\pm\pi/2.$ (9)

Defining the momentum rate ratio of two winds as:

$\displaystyle\eta$ $\displaystyle=$ $\displaystyle\frac{L_{\rm sd}/c}{4\pi
p_{0}(1+G\mid\cos\theta\mid^{m})},$ (10)

the distances of the stagnation point at the shock from the pulsar and the
star are:

$\displaystyle r_{\rm s}$ $\displaystyle=$ $\displaystyle
d\frac{\eta^{1/2}}{1+\eta^{1/2}},$ (11)

and

$\displaystyle R_{\rm s}$ $\displaystyle=$ $\displaystyle
d\frac{1}{1+\eta^{1/2}},$ (12)

respectively, with $d$ being the orbital separation.

The post-shock magnetic field is given by the Rankine-Hugoniot relations as
(Kennel & Coroniti 1984a,b):

$\displaystyle B$ $\displaystyle=$ $\displaystyle
3(1-4\sigma)\left[\frac{L_{\rm sd}\sigma}{r_{\rm
s}^{2}c(1+\sigma)}\right]^{1/2},$ (13)

and the corresponding magnetic energy density is

$\displaystyle u_{\rm B}$ $\displaystyle=$ $\displaystyle\frac{B^{2}}{8\pi},\
\rm{for}\ \phi\in[0,2\pi],$ (14)

where $\sigma$ is the magnetization parameters of the pulsar wind, which is
usually less than unity in the termination shock. Assuming $\sigma$ is
constant, we have $B\propto r_{\rm s}^{-1}$ and $u_{\rm B}\propto r_{\rm
s}^{-2}$. As the pulsar crosses the disc, the additional disc pressure will
push the shock closer to the pulsar, and therefore the post-shock magnetic
field and synchrotron radiation will be enhanced. In the following, we simply
adopt $\sigma=0.001$ throughout the text.

It is usually believed that the main target photons for IC scattering in
$\gamma$-ray binaries are provided by the luminous massive star. The energy
density of the stellar photon field at the shock is

$\displaystyle u_{\rm star}$ $\displaystyle=$ $\displaystyle\frac{L_{\rm
star}}{4\pi R_{\rm s}^{2}c},\ \rm{for}\ \phi\notin[\phi_{\rm
d,\pm}\pm\Delta\phi_{\rm d}],$ (15)

with $L_{\rm star}$ being the stellar luminosity. The contribution of stellar
photons to IC scattering is expected to peak around super-conjunction (SUPC)
or periastron phases, where the scattering is most efficient or the photon
field is densest. Alternatively, the IR emission generated by the disc may
also contribute to IC, but its contribution to IC is usually much less than
that of stellar photons (van Soelen & Meintjes 2011). Therefore, an additional
seed photon component is required to explain the two-peak profiles of TeV
flux.

Khangulyan et al. (2012) proposed that the shock heating of the disc can
supply extra soft photons for the IC process and therefore increase the
$\gamma$-ray luminosity. The amount of energy available for the seed photons
generated by the shock heating can be estimated from the kinetic energy of the
disc (Zabalza et al. 2011b; van Soelen et al. 2012). Under the dynamic
balance, the energy density of the radiation field due to the shock-heating
process may be estimated with (Chen et al. 2019)

$\displaystyle u_{\rm disc}$ $\displaystyle\simeq$ $\displaystyle\frac{\xi
L_{\rm sd}}{4\pi r_{\rm s}^{2}c},\ \rm{for}\ \phi\in[\phi_{\rm
d,\pm}\pm\Delta\phi_{\rm d}],$ (16)

where $\xi$ is the heating efficiency. In the following, we simply adopt
$\xi=0.5$ throughout the text. Its maximal contribution to IC in the shock
approximately coincides with the pulsar passing through the midplane of the
disc.

###  2.3 Radiation processes

The correlations as seen in the keV and TeV flux of $\gamma$-ray binaries
indicate a common particle population and emitting region radiating at both
energy bands (Zabalza et al. 2011a). In the pulsar scenario, the termination
shock compresses the magnetic field and accelerates electrons carried by the
pulsar wind. Shock-accelerated electrons emit synchrotron radiation in X-rays
and upscatter the soft photons to $\gamma$-rays. In the following, we adopt a
simple radiation model to minimise model parameters and simplify calculations
while maintaining the ability to reproduce the observed orbital modulations of
keV and TeV fluxes.

Similar to Dubus et al. (2017), we assume that the emitting electrons follow a
mono-energetic distribution with a Lorentz factor of $\gamma=10^{6}$, given
the fact that electrons of this energy dominate the keV and TeV emissions. We
should note that the realistic particle distributions in the shock could be
more complicated, and in a more sophisticated model, the acceleration and
cooling processes of shocked electrons should be considered. As we focus on
the orbital modulations of keV and TeV flux, assuming a more complicated
distribution of particles will not significantly change our results (e.g.
Dubus et al. 2017; Chen et al. 2019). Under the mono-energetic distribution,
the number of emitting particles is related to the particle luminosity $L_{\rm
p}$ as (Dubus et al. 2017):

$\displaystyle N_{\rm e}$ $\displaystyle=$ $\displaystyle\frac{L_{\rm
p}}{\gamma m_{\rm e}c^{2}}\times\tau_{\rm esc},$ (17)

where $m_{\rm e}$ is the electron mass, and $\tau_{\rm esc}\sim d/c$ is the
escape timescale of electrons from the emitting region. We assume that only a
fraction of spin-down luminosity goes into particle luminosity (i.e. $L_{\rm
p}/L_{\rm sd}<1$). Also, we ignore the IC cooling, which will be a good
approximation for $\gamma=10^{6}$ in Be/$\gamma$-ray binaries with larger
binary separations. The effect of radiative cooling will be important for
those binaries harbouring an O star with a much more compact orbit.

The synchrotron power radiated by an electron of Lorentz factor $\gamma$ is

$\displaystyle P_{\rm{syn}}$ $\displaystyle=$
$\displaystyle\frac{4}{3}c\sigma_{\rm T}u_{\rm B}\gamma^{2}\beta^{2},$ (18)

with the characteristic energy of

$\displaystyle\epsilon_{\rm{syn}}$ $\displaystyle=$ $\displaystyle
h\frac{q_{\rm e}B}{2\pi m_{\rm e}c}\gamma^{2}\sim 1.15B_{-1}\gamma_{6}^{2}\
{\rm keV},$ (19)

where $h$ is Plank constant, $\sigma_{\rm T}$ is the Thompson cross-section,
$q_{\rm e}$ is the charge of the electron, $\beta=\sqrt{1-\gamma^{-2}}$, and
$Q_{x}=Q/10^{x}$ is adopted.

When the pulsar is out of the disc region (i.e. $\phi\notin[\phi_{\rm
d,\pm}\pm\Delta\phi_{\rm d}]$), the seed photons for IC scattering mainly come
from the black-body photons emitted by the star. As stellar photons are
radiated radially, anisotropic scattering should be considered. The
anisotropic IC scattering power for a single electron is (Dubus & Cerutti
2013)

$\displaystyle P_{\rm{IC}}$ $\displaystyle\simeq$
$\displaystyle\frac{4}{3}c\sigma_{\rm
T}u_{\rm{star}}\left(1-\beta\mu\right)\left[\left(1-\beta\mu\right)\gamma^{2}-1\right]F_{\rm{KN}}\left(\gamma\right),$
(20)

where $\mu=\@vec{e}_{\rm{obs}}\cdot\@vec{e}_{\rm{psr}}$ is the scattering
angle, and $F_{\rm KN}(\gamma)$ is the reduction factor due to Klein-Nishina
effect (Moderski et al. 2005; Khangulyan et al. 2014). During the disc
passages (i.e. $\phi\in[\phi_{\rm d,\pm}\pm\Delta\phi_{\rm d}]$), seed photons
are generated by shock heating of the disc. In this case, we assume that the
scattering process is isotropic, and therefore

$\displaystyle P_{\rm{IC}}$ $\displaystyle\simeq$
$\displaystyle\frac{4}{3}c\sigma_{\rm
T}u_{\rm{disc}}\gamma^{2}\beta^{2}F_{\rm{KN}}\left(\gamma\right).$ (21)

The characteristic energies of up-scattered photons can be estimated as

$\displaystyle\epsilon_{\rm IC}$ $\displaystyle\simeq$
$\displaystyle\left\\{\begin{array}[]{ll}4\gamma^{2}\epsilon_{0}\sim
4.0\gamma_{6}^{2}\epsilon_{0,\rm eV}\ {\rm TeV},&\hbox{for
$\gamma\epsilon_{0}/m_{\rm{e}}c^{2}\ll 1$,}\\\ \gamma m_{\rm{e}}c^{2}\sim
0.5\gamma_{6}\ {\rm TeV},&\hbox{for $\gamma\epsilon_{0}/m_{\rm{e}}c^{2}\gg
1$,}\end{array}\right.$ (24)

where $\epsilon_{0}$ is the characteristic energy of seed photons.

The keV X-ray luminosity due to synchrotron radiation is given by

$\displaystyle L_{\rm X}$ $\displaystyle\sim$ $\displaystyle N_{\rm e}\times
P_{\rm syn},$ (25)

and the TeV $\gamma$-ray luminosity due to IC scattering is given by

$\displaystyle L_{\rm\gamma}$ $\displaystyle\sim$ $\displaystyle N_{\rm
e}\times P_{\rm IC}.$ (26)

With a constant injection rate of emitting particles, the modulations of
synchrotron and IC scattering luminosities are mainly caused by variations in
the energy densities of the magnetic field and photon field at the shock. In
addition, the anisotropy of IC will also slightly affect the $\gamma$-ray
flux.

##  3 Results

In this section, we compare our calculated light curves with the observational
data of four detected Be/$\gamma$-ray binaries. In the model, the keV and TeV
flux are produced by synchrotron radiation and IC scattering in the shock,
respectively, as depicted above.

###  3.1 PSR B1259-63/LS 2883

Figure 3:  Orbital geometry of PSR B1259-63/LS 2883. The black pentagrams are
the phases of periastron and apastron, and superior and inferior conjunctions,
while the empty circles mark the time intervals of every 20 d from periastron.
The positions of double peaks in keV/TeV light curves are shown in green. The
inclined disc projected on the orbital plane is illustrated in orange, with
the dashed line being the midplane of the disc.

Figure 4:  Calculated X-ray (top) and $\gamma$-ray (bottom) light curves of
PSR B1259-63/LS 2883 with comparisons of observational data. The X-ray data
during the 2004, 2007, 2010, 2014, and 2017 periastron passages are taken from
Chernyakova et al. (2006, 2009, 2014, 2015) and Tam et al. (2018),
respectively. The H.E.S.S. data are from Romoli et al. (2017).

PSR B1259-63/LS 2883 is the first $\gamma$-ray binary where the compact object
was identified as a rotational pulsar, and it is also the first variable
Galactic source detected in TeV energies (Johnston et al. 1992; Aharonian et
al. 2005). The massive companion, LS 2883, with a bolometric luminosity of
$L_{\star}\simeq 2.2\times 10^{38}\ {\rm erg\ s^{-1}}$, is characterised by an
equatorial disc inclined to the orbital plane. Long-term radio timing
observations of PSR B1259-63 allow for precise measurements of related binary
parameters, including the orbital period $P_{\rm o}=1236.724526\ \rm{d}$, the
eccentricity $e=0.8698797$, the inclination angle $i_{\rm o}=154^{\circ}$ ,
and the argument of periastron $\omega_{\rm p}=138^{\circ}_{.}665013$ (Wang et
al. 2004; Shannon et al. 2014; Miller-Jones et al. 2018). According to the
parallax data of the Gaia DR2 Archive, the distance was updated to be $d_{\rm
L}=2.39\pm 0.19\ {\rm kpc}$ (Gaia Collaboration 2018; Chernyakova et al.
2020a). The orbit of PSR B1259-63/LS 2883 is illustrated in Fig. 3.

Since its discovery in 1992, the binary has been intensively monitored by many
telescopes, which collect a huge amount of data. Several simultaneous multi-
wavelength campaigns have been performed on the system (Chernyakova et al.
2014, 2015, 2020a, 2021). One of the most noticeable features of PSR
B1259-63/LS 2883 is the two asymmetrical peaks as seen in its keV and TeV
light curves. Specifically, the X-ray flux shows a rapid increase around
$T_{0}-15\ {\rm d}$, and reaches its maximum about $5\ {\rm d}$ later (where
$T_{0}$ corresponds to the periastron epoch). After that, the flux decays
slowly and is then followed by another peak around $T_{0}+20\ {\rm d}$ (e.g.
Tam et al. 2018). The $\gamma$-ray flux detected by the HESS telescope
combined with all available data during the past decade also displays a
similar correlated modulation as seen in the X-ray band (H.E.S.S.
Collaboration 2020). The phase ranges of two peaks as seen in keV and TeV
light curves are shown in Fig. 3 in green.

The bimodal structures of keV/TeV flux are thought to be caused by the pulsar
wind interacting with stellar outflow, in particular with the inclined disc
(Chen et al. 2019, and references therein). In Fig. 4, we use the emission
model described in previous sections to calculate the light curves and compare
them with observational data. The momentum rate of the polar wind adopted in
the calculation is $p_{0}=6.3\times 10^{26}\ {\rm g\ cm\ s^{-2}}$,
corresponding to a typical mass-loss rate of $\dot{M}\sim 5\times
10^{-8}M_{\odot}\ {\rm yr}^{-1}$ and velocity of $v_{\rm w}\sim 2\times
10^{8}\ {\rm cm\ s^{-1}}$ for the polar wind of Be stars. The disc-to-wind-
pressure contrast is $G=50$, which is estimated from the observed amplitudes
of the peaks. As the pulsar moves into the disc region, the size of the shock
shrinks because of the additional disc pressure. As $u_{\rm B}\propto r_{\rm
s}^{-2}$ and $u_{\rm disc}\propto r_{\rm s}^{-2}$, the energy densities of
magnetic field and photon field are increased during disc passages. Therefore,
this process will enhance synchrotron radiation and IC scattering, producing
two peaks as seen in the keV and TeV light curves. The second peak is higher
than the one before periastron because of the Doppler boosting effect when the
pulsar moves around INFC. The shocked flow velocity adopted in calculations is
$v_{\rm s}=c/3$ (Dubus et al. 2010). The disc parameters can be estimated by
the fitting the observational data: (1) the positions of the two peaks suggest
that the midplane of the disc is located at $\phi_{\rm d,\pm}=12^{\circ}\pm
90^{\circ}$; (2) the widths of the peaks indicate that the disc opening angle
projected on the orbital plane is about $\Delta\phi_{\rm d}\simeq
11.8^{\circ}$, with an inclination angle of $i_{\rm d}\simeq 145^{\circ}$.
This agrees with the measurement by Shannon et al. (2014), who find that the
angle between the stellar axis of LS 2883 and the orbital axis is about
$35^{\circ}$ (we note that PSR B1259-63 moves clockwise on its orbit as shown
by Miller-Jones et al. 2018). The position of the disc projected on the
orbital plane is also illustrated in Fig. 3.

We note that the observed X-ray flux decays somehow slower than the one
predicted by the model. Furthermore, Chernyakova et al. (2021) discovered a
third X-ray peak during the 2021 periastron passage which has never been
detected before. The peak started to rise around 30 d after periastron, and
did not show any clear correlated emission at other wavelengths. The slow
decay behaviour and the additional peak could be due to the disc matter being
piled up at the shock after disc passages, as discussed in Chen et al. (2019).
Alternatively, the inhomogeneity and the clumps in the stellar outflows can
also lead to the enhancement of X-rays (Chernyakova et al. 2021).

###  3.2 PSR J2032+4127/MT91 213

Table 1: Three orbital solutions of PSR J2032+4127/MT91 213 given by Ho et al. (2017).  Parameters | Model 1 | Model 2 | Model 3
---|---|---|---
$e$ | 0.936 | 0.961 | 0.989
$P_{\rm o}(\rm{d})$ | 16000 | 17000 | 17670
$T_{0}$(MJD) | 58053 | 58069 | 58068
$\omega_{\rm p}(^{\circ})$ | 52 | 40 | 21
$a\sin{i}$ (lt-s) | 7138 | 9022 | 16335

Figure 5:  Orbital geometry of PSR J2032+4127/MT91 213 around periastron with
the orbital solution of model 2 in Table 1.

(a) Model 1

(b) Model 2

(c) Model 3

Figure 6: Calculated X-ray and $\gamma$-ray light curves of PSR
J2032+4127/MT91 213 under three different orbital solutions of Ho et al.
(2017). The X-ray data are analysed in this work, while the $\gamma$-ray data
are taken from Abeysekara et al. (2018). The inserted plots display the light
curves from between $T_{0}-120$ d to $T_{0}+120$ d.

PSR J2032+4127 was firstly identified as an isolated radio-loud $\gamma$-ray
pulsar with a spin period of $P=143.2\ {\rm ms}$ (Abdo et al. 2009; Camilo et
al. 2009). Later radio observations found that the pulsar displays an
extraordinary increase in the spin-down rate, which was attributed to the
pulsar being in orbit with a massive star (Lyne et al. 2015). The multi-
wavelength monitoring conducted by Ho et al. (2017) further confirmed that PSR
J2032+4127 is moving in a highly eccentric orbit ($e\sim 0.94-0.99$) around
the B0Ve star MT91 213 with a very long period $P_{\rm o}\sim 16000-17670\
{\rm d}$. The above characteristics suggest that PSR J2032+4127/MT91 213 is a
similar $\gamma$-ray binary to PSR B1259-63/LS 2883.

The X-ray flux of PSR J2032+4127/MT91 213 increased steadily before October
2017, and was then accompanied by a rapid dip around periastron (Li et al.
2017, 2018). After that, the flux was followed by another flare that lasted
for several tens of days (Coe et al. 2019; Pal et al. 2019; Ng et al. 2019).
The VERITAS and MAGIC telescopes also detected TeV $\gamma$-rays from the
binary (Abeysekara et al. 2018). The TeV flux was found to increase steadily
before October 2017, and was then followed by a short dip at periastron. A few
days later, the flux increased again with a level comparable to the pre-
periastron phase.

Due to its long period and high eccentricity, the orbital parameters of PSR
J2032+4127/MT91 213 are not well measured. Ho et al. (2017) presented three
different binary models, and the related orbital parameters are summarised in
Table 1. The orbit of model 2 is illustrated in Fig. 5.

Given the similar nature of the two sources, we expect the emissions from PSR
J2032+4127/MT91 213 can be explained by a similar model as that explaining PSR
B1259-63/LS 2883. Therefore, we apply the inclined disc model to calculate the
expected orbital modulations of keV and TeV flux with the orbital solutions
provided by Ho et al. (2017). The comparisons with observational data around
periastron are presented in Fig. 6. The X-ray data are reproduced from the
Swift/XRT website222https://www.swift.ac.uk/user_objects/ covering the time
from MJD 54087 to MJD 58560, while the $\gamma$-ray data are taken from
Abeysekara et al. (2018). The model parameters adopted in calculations are the
same as PSR B1259-63/LS 2883, except for the inclination and position angle of
the disc (i.e. $i_{\rm d}$ and $\phi_{\rm d}$), and the shocked flow velocity
adopted here is $v_{\rm s}\sim 0.1c$. As we can see, although the predicted
light curves following model 1 agree well with the data far before periastron,
the calculated flux around periastron does not match the observations. As for
model 3 with a very high eccentricity ($e=0.989$), the model also fails to fit
the overall observational data. Instead, our calculation results suggest that
model 2 of Ho et al. (2017) with modest values of eccentricity and period
($e=0.961,P_{\rm o}\simeq 17000\ {\rm d}$) provides a reasonable fit to the
data, especially around the periastron passage. Our fitting results suggest
that the position of the midplane of the disc is located at $\phi_{\rm
d,\pm}=165^{\circ}\pm 90^{\circ}$ with a rather small inclination angle
$i_{\rm d}\simeq 9^{\circ}$. This corresponds to a relatively large disc
opening angle on the orbital plane with $\Delta\phi_{\rm d}\simeq
48.6^{\circ}$. The position of the disc is displayed in Fig. 5.

The observed X-ray flux far before periastron is slightly higher than the
predicted light curve of model 2. This is probably due to oversimplifications
in our disc model. It has been suggested that the disc opening angle becomes
larger with growing radius (Carciofi & Bjorkman 2006, 2008), and therefore the
pulsar would interact with the disc far before periastron. This explains why
the observed X-ray flux is higher than the model prediction. Alternatively,
the discrepancy could also be due to the inhomogeneity of stellar outflow or
an additional wind driven by the disc.

###  3.3 HESS J0632+057

HESS J0632+057 was discovered by the survey around the Monoceros region
(Aharonian et al. 2007). Its optical counterpart, MWC 148, is a B0pe-type star
with a bolometric luminosity of $L_{\star}\simeq 1.6\times 10^{38}\ {\rm erg\
s^{-1}}$. Although there is no direct detection of pulsed signals, the compact
object of HESS J0632+057 has been widely believed to be a rotational NS (Yi &
Cheng 2017; Bosch-Ramon et al. 2017; Barkov & Bosch-Ramon 2018; Malyshev et
al. 2019).

The system shows similar spectral index and flux variability as seen in other
Be/$\gamma$-ray binaries. In particular, the keV/TeV light curves are
characterised by two asymmetric peaks. The primary peak around phase $0.3-0.4$
shows a sharp increase and decrease, while the second one around $0.6-0.9$ has
a flat shape with a longer duration. A clear dip around phase $0.4-0.5$ was
seen at both keV and TeV bands, and the X-ray light curves also exhibit a
plateau around $0.1-0.2$ (Bongiorno et al. 2011; Aliu et al. 2014; Malyshev et
al. 2019; Adams et al. 2021; Tokayer et al. 2021).

Unfortunately, the orbital parameters of HESS J0632+057 are not yet well
measured, and different studies give different results. Based on optical
spectroscopy of the Be companion MWC 148, Casares et al. (2012) proposed an
orbital period of $P_{\rm o}=321\pm 5\ \rm{d}$ with an eccentricity of
$e=0.83\pm 0.08$ and a periastron phase of $\Phi_{\rm p}=0.967\pm 0.008$.
Alternatively, a slightly shorter period of $P_{\rm o}=313^{+11}_{-9}\ \rm{d}$
with a less elliptical orbit ($e=0.643\pm 0.29$) was derived by Moritani et
al. (2018) according to modulation of the H$\alpha$ emission line. It was also
suggested that the periastron phase is around $\Phi_{\rm p}=0.663$, which is
located on the opposite side of the orbit suggested by Casares et al. (2012).
Recently, Malyshev et al. (2019) and Adams et al. (2021) refined the orbital
period to $317.3\pm 0.7\ \rm{d}$ using the long-term X-ray observational data
of Swift/XRT, XMM-Newton, Chandra, NuSTAR, and Suzaku. Based on the positions
and the relative widths of the X-ray peaks, Malyshev et al. (2019) suggested
that the periastron is located around $\Phi_{\rm p}\sim 0.4$, with an orbital
eccentricity of $e\sim 0.5$. In Table 2, we summarise the orbital solutions of
Casares et al. (2012) and Moritani et al. (2018), and the corresponding orbits
are illustrated in Fig. 7. The phase ranges of the two peaks as seen in keV
and TeV light curves are also marked in green.

The similar modulations of keV and TeV flux between HESS J0632+057 and PSR
B1259-63/LS 2883 indicate that HESS J0632+057 may also host a non-accreting
NS, and the two-peak profiles are also caused by interactions between the
presumptive pulsar and the inclined disc. Therefore, in Fig. 8, we calculate
the expected shock emissions under the orbital models as given in Table 2, and
compare them with the observational data. The X-ray data are taken from
Malyshev et al. (2019), while the $\gamma$-ray data are adopted from Aliu et
al. (2014). The observational data are refolded with different orbital periods
in Table 2 as

$\displaystyle\Phi$ $\displaystyle=$
$\displaystyle\frac{T-T_{0}}{P_{\mathrm{o}}}-\mathrm{int}\left(\frac{T-T_{0}}{P_{\mathrm{o}}}\right),$
(27)

where $T_{0}={\rm MJD}\ 54857.0$ is the zero phase time (Bongiorno et al.
2011). As displayed in Fig. 8, the expected keV and TeV emissions under the
inclined disc model with the orbit solutions of Casares et al. (2012) and
Moritani et al. (2018) fail to reproduce the observed modulations. This can
also be seen from the orbital phase ranges of two peaks as illustrated in the
upper and middle panels of Fig. 7. If the double humps are indeed caused by
the presumptive pulsar passing through the disc, the true anomaly interval
between two peaks should be $180^{\circ}$, and therefore even by changing the
values of $i_{\rm d}$ and $\phi_{\rm d}$, it is still difficult to match the
data.

Here, we propose an alternative orbital geometry for HESS J0632+057. Because
the pulsar moves faster at periastron and slower at apastron, the periastron
phase should be located at the shortest separation between the peaks, around
phases $0.3-0.7$ (Malyshev et al. 2019). Also, the primary peak is narrower
than the second one, which puts further constraints on the periastron phase
around $0.3-0.4$. Therefore, we assume that the periastron phase is around
$\Phi_{\rm p}\sim 0.35$. Also, the positions and separation of the two peaks
suggest a less eccentric orbit with $e\sim 0.35$. Although this eccentricity
is much smaller than that of Casares et al. (2012), it is still within the
error range of Moritani et al. (2018). The suggested orbital parameters are
summarised as model 3 of Table 2, and the corresponding orbit and expected
light curves are shown in the bottom panels of Figs. 7 and 8, respectively.
The overall double-peak profiles in the keV and TeV light curves are
relatively consistent with calculated results under the orbital solution of
model 3, except for the X-ray plateau before the primary peak. The primary
peak is significantly higher than the second one because the pulsar crosses
the inner part of the disc near periastron, which pushes the shock closer to
the pulsar. As $u_{\rm B}\propto r_{\rm s}^{-2}$ and $u_{\rm disc}\propto
r_{\rm s}^{-2}$, enhanced synchrotron radiation and IC scattering are expected
especially during the first disc passages. Our fitting results suggest that
the position of the midplane and the inclination angle of the disc are
$\phi_{\rm d,\pm}=72^{\circ}\pm 90^{\circ}$ and $i_{\rm d}\simeq 12^{\circ}$,
respectively. The corresponding disc opening angle on the orbital plane is
$\Delta\phi_{\rm d}\simeq 34.5^{\circ}$. The derived values of the orbital
parameters and the position of the disc are roughly consistent with the
analysis of Malyshev et al. (2019).

It is necessary to point out that our calculated X-ray light curve does not
closely match the X-ray plateau. The plateau could be caused by the clumps of
the stellar outflows, as in the case of PSR B1259-63, because the clumpy
outflows can push the shock closer to the pulsar side, causing a stronger
magnetic field and therefore the enhanced X-ray emission (Chernyakova et al.
2021). Alternatively, the plateau could be due to the Doppler-boosting effect
if the inferior conjunction is located around phase 0.1-0.2. However, the lack
of a plateau phase in the TeV light curve makes this scenario unlikely.

###  3.4 LS I +61∘303

LS I +61∘303 was firstly discovered by the Cos B satellite as an unknown
bright $\gamma$-ray emitter (Hermsen et al. 1977), and later Gregory & Taylor
(1978) reported variable radio emissions from the exact location. Optical
spectroscopy suggested the system consists of a stellar-mass compact object in
orbit with an early B0Ve-type massive star (Paredes & Figueras 1986; Casares
et al. 2005). The nature of the compact object remains unclear even after more
than 40 years of multi-wavelength observations (e.g. see Marcote 2017 and
references therein). For quite a long time, LS I +61∘303 has been considered
as a micro-quasar (e.g. Massi, Ros & Zimmermann 2012; Massi & Torricelli-
Ciamponi 2014; Jaron 2021); however, the lack of an accretion signal in X-rays
suggests that the non-thermal emission is more likely powered by spin-down of
a young pulsar (Maraschi & Treves 1981; Harrison et al. 2000; Leahy 2004;
Zdziarski et al. 2010; Zabalza et al. 2011a). Although several campaigns were
conducted to search for pulsed radio and X-ray signals, no pulsations were
found (Sidoli et al. 2006; McSwain et al. 2011; Cañellas et al. 2012).
Recently, Weng et al. (2021) reported the detection of a transient periodic
signal of 269.196 ms from LS I +61∘303 on 2020 Jan 7 with the FAST telescope.
However, no more pulsations were found in other FAST observation campaigns.
This is probably caused by the strong free-free absorption by the stellar
outflows, in particular with the presence of the stellar disc (Chen et al.
2021a).

According to a Bayesian analysis of the radio data, the orbital period of LS I
+61∘303 is $26.4960\pm 0.0028$ d (Gregory 2002). The most peculiar feature of
this system is its superorbital modulations with $P_{\rm sup}\simeq 1628\pm
48$ d as seen at all wavelengths, which has not been found for other
$\gamma$-ray binaries (e.g., Li et al. 2011; Chernyakova et al. 2012, 2017;
Ackermann et al. 2013; Ahnen et al. 2016; Massi & Torricelli-Ciamponi 2016;
Jaron et al. 2018). The origin of its long-term modulations remains unknown,
although several models have been proposed, including precession of the Be
disc or the relativistic jet (Massi & Torricelli-Ciamponi 2014; Saha et al.
2016; Xing et al. 2017; Jaron 2021).

A simultaneous observational campaign with Swift/XRT, XMM-Newton, and the
MAGIC telescope was conducted by Anderhub et al. (2009), which revealed a
significant correlation between the keV and TeV flux. The prominent peaks
around phase 0.65 are clearly seen at both energy bands, where the orbital
phase is defined as Eq. (27) with $T_{0}={\rm JD}\ 2443366.775$ (Gregory
2002). Following the short dip around $0.7-0.8$, the X-ray flux displayed
another slow increase from phase 0.8 to 1.0. This similar behaviour is also
being found in the TeV band.

Unfortunately, similar to HESS J0632+057, the binary parameters of LS I
+61∘303 are not well measured, and several orbital solutions are derived by
different methods (Hutchings & Crampton 1981; Casares et al. 2005; Grundstrom
et al. 2007; Aragona et al. 2009; Kravtsov et al. 2020). The optical
spectroscopy of Casares et al. (2005) suggested the eccentricity is $e=0.72\pm
0.15$ with periastron occurring around phase $\Phi_{\rm p}=0.23\pm 0.02$,
while the H$\alpha$ monitoring campaign of Grundstrom et al. (2007) suggested
$e=0.55\pm 0.05$ with $\Phi_{\rm p}=0.24\pm 0.04$. More recently, a less
elliptical orbit ($e<0.2$) with periastron around $\Phi_{\rm p}\sim 0.6$ was
derived by Kravtsov et al. (2020) based on optical linear polarisation
variability curves. In Table 3, we summarise the representative orbital models
of Casares et al. (2005) and Grundstrom et al. (2007), and corresponding
orbits are illustrated in Fig. 9. The phase ranges of the two peaks as seen in
keV and TeV bands are also marked in green.

In Fig. 10, we present the corresponding expected keV and TeV light curves
under the inclined disc model with the orbital models of Table 3. As LS I
+61∘303 shows long-term modulations in its multi-wavelength emissions, we
adopt the simultaneous X-ray and $\gamma$-ray data from Anderhub et al.
(2009). As presented in the top and middle panels of Fig. 10, the predicted
light curves under the orbit solution of Casares et al. (2005) and Grundstrom
et al. (2007) fail to reproduce the observed modulations of keV/TeV flux. This
can also be seen from the orbital phase ranges of the two peaks as shown in
Fig. 9. Similar to the analysis in HESS J0632+057, if the double-hump
structures are indeed due to the pulsar–disc interactions, then the periastron
phase should be located at the shortest separation between the peaks around
phase $0.7-0.8$. Also, the position of the prominent peak further indicates
that the periastron phase would be around $0.6-0.7$. Therefore, we assume that
$\Phi_{\rm p}\sim 0.7$. This value is roughly consistent with that of Kravtsov
et al. (2020). The orbital eccentricity with $e\sim 0.40$ is derived from the
relative width of the two peaks, and is slightly higher than that of Kravtsov
et al. (2020). The suggested orbital parameters are summarised as model 3 of
Table 3, and the corresponding orbit and expected keV/TeV light curves are
shown in the bottom panels of Figs. 9 and 10, respectively. Our result
suggests that the position of the midplane and the inclination angle of the
disc are $\phi_{\rm d,\pm}=30^{\circ}\pm 90^{\circ}$ and $i_{\rm d}\simeq
10^{\circ}$, respectively. These correspond to the disc opening angle on the
orbital plane with $\Delta\phi_{\rm d}\simeq 42.5^{\circ}$. Similar to HESS
J0632+057, LS I +61∘303 also exhibits an X-ray plateau around phase 0.45-0.55,
which could be caused by the inhomogeneity of the stellar outflows.

Table 2: Orbital solutions of HESS J0632+057.

Parameters | Model 1† | Model 2‡ | Model 3§
---|---|---|---
$e$ | $0.83\pm 0.08$ | $0.643\pm 0.29$ | 0.35
$P_{\rm o}(\rm{d})$ | $321\pm 5$ | $313^{+11}_{-8}$ | $317.3$
$\Phi_{\rm p}$ | $0.967$ | $0.663$ | $0.35$
$\omega_{\rm p}(^{\circ})$ | $129\pm 17$ | $271\pm 29$ | $129$

${\dagger}$ Casares et al. 2012;
${\ddagger}$ Moritani et al. 2018;
$\S$ This work. The values of the orbital period $P_{\rm o}$ and longitude of
periastron $\omega_{\rm p}$ are taken from Adams et al. (2021) and Casares et
al. (2009), respectively. The eccentricity and periastron phase adopted here
are slightly different from that of Malyshev et al. (2019).

(a) Model 1

(b) Model 2

(c) Model 3

Figure 7: Orbital geometries of HESS J0632+057 with three different orbital
solutions given in Table 2.

(a) Model 1

(b) Model 2

(c) Model 3

Figure 8: Calculated light curves of HESS J0632+057 with comparisons of
observational data under three different orbital models as given in Table 2.
The X-ray data are taken from Malyshev et al. (2019), and the $\gamma$-ray
data are taken from Aliu et al. (2014). We note that the data are refolded
with different orbital periods via Eq. (27).

Table 3: Orbital solutions of LS I +61∘303\.

Parameters | Model 1${\dagger}$ | Model 2${\ddagger}$ | Model 3$\S$
---|---|---|---
$e$ | $0.72\pm 0.15$ | $0.55\pm 0.05$ | $0.40$
$P_{\rm o}(\rm{d})$ | $26.496$ | $26.496$ | $26.496$
$\Phi_{\rm p}$ | $0.23\pm 0.02$ | $0.301\pm 0.011$ | 0.70
$\omega_{\rm p}(^{\circ})$ | $21\pm 13$ | $57\pm 9$ | $40.5$

${\dagger}$ Casares et al. 2005;
${\ddagger}$ Grundstrom et al. 2007;
$\S$ This work. The value of longitude of periastron $\omega_{\rm p}$ is taken
from Aragona et al. 2009.

(a) Model 1

(b) Model 2

(c) Model 3

Figure 9: Orbital geometries of LS I +61∘303 with three different orbital
solutions given in Table 3.

(a) Model 1

(b) Model 2

(c) Model 3

Figure 10: Calculated light curves of LS I +61∘303 with comparisons of
observational data under three different orbital models as given in Table 3.
The X-ray and $\gamma$-ray data are taken from Anderhub et al. (2009).

Table 4: Parameters of PSR B1259-63/LS 2883, PSR J2032+4127/MT91 213, HESS
J0632+057, and LS I +61∘303\.

Parameters | PSR B1259-63a | PSR J2032+4127b | HESS J0632+057c | LS I +61∘303d
---|---|---|---|---
observational parameters |  |  |  |
$e$ | 0.8698797 | 0.961 | 0.35${\dagger}$ | 0.40${\dagger}$
$\Phi_{\rm p}$ | - | - | 0.35${\dagger}$ | 0.70${\dagger}$
$P(\rm{d})$ | 1236.724526 | 17000 | 317.3 | 26.496
$d_{\rm L}(\rm{kpc})$ | 2.39 | 1.33 | 1.60 | 2.30
$i_{\rm o}(^{\circ})$ | 154 | $\sim 60$ | $\sim 60$ | $\sim 45$
$\phi_{\rm o}(^{\circ})$ | 132 | $\sim 229$ | $\sim 321$ | $\sim 50$
$L_{\rm sd}(\rm erg\ s^{-1})$ | $8.2\times 10^{35}$ | $1.7\times 10^{35}$ | $5\times 10^{35}{\dagger}$ | $5\times 10^{35}{\dagger}$
$L_{\rm star}(\rm erg\ s^{-1})$ | $2.2\times 10^{38}$ | $2.6\times 10^{38}$ | $1.6\times 10^{38}$ | $2.6\times 10^{38}$
model parameters |  |  |  |
$p_{0}(\rm{g\cdot cm\cdot s^{-2}})$ | $6.3\times 10^{26}$ | $6.3\times 10^{26}$ | $2.5\times 10^{26}$ | $1.3\times 10^{26}$
$G$ | $50$ | 50 | 50 | 50
$m$ | $100$ | 100 | 100 | 100
$i_{\rm d}(^{\circ})$ | $145$ | $9$ | $12$ | $10$
$\phi_{\rm d}(^{\circ})$ | $12$ | $165$ | $72$ | $30$
$\Delta\phi_{\rm d}(^{\circ})$ | $11.8$ | $48.6$ | $34.5$ | $42.5$
$\sigma$ | $10^{-3}$ | $10^{-3}$ | $10^{-3}$ | $10^{-3}$
$\xi$ | 0.5 | 0.5 | 0.5 | 0.5
$v_{\rm s}(c)$ | 1/3 | 0.1 | 0.02 | 0.1

References:
$a$ Negueruela et al. 2011; Shannon et al. 2014; Miller-Jones et al. 2018.
$b$ Camilo et al. 2009; Lyne et al. 2015; Ho et al. 2017.
$c$ Casares et al. 2012; Moritani et al. 2018; Malyshev et al. 2019; Adams et
al. 2021.
$d$ Gregory 2002; Aragona et al. 2009.
${\dagger}$ This work.

##  4 Discussion and Conclusion

Be/$\gamma$-ray binaries provide information as to the relationship between Be
stars and high-energy astrophysics (Lamberts 2016; Moritani & Kawachi 2021).
Interaction between these compact objects and their stellar outflow is
believed to play an essential role in shaping their multi-wavelength light
curves. In this paper, we present our study of the correlation between the keV
and TeV emission of Be/$\gamma$-ray binaries under the pulsar scenario. We
show that the double-hump structures of keV/TeV flux are caused by the pulsar
wind–stellar disc interactions, during which synchrotron and IC radiations are
enhanced due to the increases in magnetic field and photon field in the shock.
Within this scenario, the disc parameters can be estimated by fitting the
light curves.

We found that the keV/TeV light curves of PSR B1259-63/LS 2883 and PSR
J2032+4127/MT91 213 can be well explained by the inclined disc model. Our
fitting results suggest that the disc midplane is located at $\phi_{\rm
d,\pm}=12^{\circ}\pm 90^{\circ}$ ($\phi_{\rm d,\pm}=165^{\circ}\pm
90^{\circ}$) with an inclination angle of $i_{\rm d}\simeq 145^{\circ}$
($i_{\rm d}\simeq 9^{\circ}$) for LS 2883 (MT91 213). As for HESS J0632+057
and LS I +61∘303, modelling their orbital modulations with the current
available orbital solutions is challenging under the inclined disc model. We
used the shapes of light curves and the separations between two peaks to
estimate the periastron phases and the orbital eccentricities, and propose
alternative orbital geometries. We suggest that the periastron phase of HESS
J0632+057 (LS I +61∘303) is located around $\Phi_{\rm p}\sim 0.35$ ($\Phi_{\rm
p}\sim 0.70$ ) with an eccentricity of $e\sim 0.35$ ($e\sim 0.40$). We derived
the positions and inclination angles of the discs by fitting the keV/TeV light
curves (i.e. $i_{\rm d}=12^{\circ},\phi_{\rm d}=72^{\circ}$ and $i_{\rm
d}=10^{\circ},\phi_{\rm d}=30^{\circ}$ for HESS J0632+057 and LS I +61∘303,
respectively). The related observational and model parameters of four
Be/$\gamma$-ray binaries are summarised in Table 4.

Although our model can roughly fit the correlated keV/TeV light curves of
observed Be/$\gamma$-ray binaries, the discrepancies between the observational
data and model predictions suggest that there are still several limitations to
our emission model. Firstly, we assume that the emitting particles have a
mono-energetic distribution, while it is generally believed that the electron
injection spectrum follows a power-law function, and the related cooling
processes will also reshape the distribution (Khangulyan et al. 2007; Zabalza
et al. 2011a; Chen et al. 2019). Therefore, in a more sophisticated model, the
related acceleration and cooling processes of emitting particles should be
taken into account. Secondly, the shock radiations are calculated under a
simple one-zone model. It should be noted that, in the realistic case, the
emitting region could be much more complicated. Generally, it is believed that
the shock has a bow-shaped geometry wrapping around the pulsar, and the
shocked flows moving from the apex to the tail could result in more
complicated boosted emissions (Chen et al. 2021b). Also, the hydrodynamic
simulations suggested that the pulsar’s orbital motion can create another
shock on the far side of the binary, which may contribute to observed
emissions, in particular for those $\gamma$-ray binaries with more compact
orbits (Bosch-Ramon et al. 2012, 2015; Huber et al. 2021a,b; Barkov & Bosch-
Ramon 2021). Thirdly, we assume that the ram pressure of Be outflow follows a
simple $(1+G\cos\theta^{m})$ distribution, and therefore that the properties
of the disc can be simply characterised by the confinement parameter $m$ and
equator-to-pole pressure contrast $G$. In the widely adopted Gaussian
distribution for the disc density, the disc opening angle becomes larger with
growing radius (e.g. Carciofi & Bjorkman 2006, 2008). Using the Gaussian disc
model to obtain the shock position requires further assumptions about the disc
speeds (including the radial, azimuthal, and vertical directions), which are
still quite uncertain for Be stars (e.g. see Sect 4.3 of Torres et al. 2012
for a discussion). Furthermore, if the disc is a Keplerian disc, the rotation
of the disc can significantly deflect the shock structure, which has not yet
been fully investigated. The rotation of the disc can significantly change the
shock structure, because the disc pressure will push the shock apex much
closer to the pulsar if the rotation of the disc is retrograde with the
pulsar’s orbit, and further away from the pulsar if the rotation of the disc
is prograde. The Keplerian disc can also significantly deflect the tail of the
shock from the orbital plane. Depending on the angle between the shock
structure and the line of sight, the Doppler boosting effect can also affect
the observed flux. In addition, the inhomogeneity of stellar outflow will
cause additional flux variabilities (Bosch-Ramon 2013; Paredes-Fortuny et al.
2015; de la Cita et al. 2017). Also, the anisotropy of the pulsar wind might
also affect the radiation (Kong et al. 2012; Bosch-Ramon 2021). We also ignore
the effect of pair absorption on TeV flux, and this process might reduce the
TeV photons, especially around periastron and SUPC (Dubus 2006b). The ensuing
cascade emissions from the secondary electron–positron pairs will also
contribute to the observed emissions at lower energies (e.g. Bosch-Ramon et
al. 2008; Takata et al. 2017).

Nevertheless, our simple emission model is still able to capture the effects
of Be discs on the shock radiations, and the overall modulations of correlated
keV/TeV flux from the four detected Be/$\gamma$-ray binaries are reproduced
under the inclined disc model. The results could be beneficial for future
measurements of orbital parameters and searches for potential radio pulsations
from compact objects.

###### Acknowledgements.

We thank Dr. Denys Malyshev for providing us the long-term X-ray data of HESS
J0632+057 and the referee for valuable comments on the manuscript. We also
want to thank Prof. Kwong-Sang Cheng and Prof. Yun-Wei Yu for useful
discussions. This work is supported by the National Key R&D Program of China
(Grant No. 2020YFC2201400), the National Natural Science Foundation of China
(Grant No. U1838102), and the China Postdoctoral Science Foundation (Grant No.
2020M682392).

## References

  * Abdo et al. (2009) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, Science, 325, 840.
  * Abeysekara et al. (2018) Abeysekara, A. U., Benbow, W., Bird, R., et al. 2018, ApJ, 867, L19.
  * Ackermann et al. (2013) Ackermann, M., Ajello, M., Ballet, J., et al. 2013, ApJ, 773, L35.
  * Adams et al. (2021) Adams, C. B., Benbow, W., Brill, A., et al. 2021, arXiv e-prints, arXiv:2109.11894.
  * Aharonian et al. (2005a) Aharonian, F., Akhperjanian, A. G., Aye, K.-M., et al. 2005a, Science, 309, 746.
  * Aharonian et al. (2005b) Aharonian, F., Akhperjanian, A. G., Aye, K.-M., et al. 2005b, A&A, 442, 1.
  * Aharonian et al. (2007) Aharonian, F. A., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2007, A&A, 469, L1.
  * Ahnen et al. (2016) Ahnen, M. L., Ansoldi, S., Antonelli, L. A., et al. 2016, A&A, 591, A76.
  * Albert et al. (2006) Albert, J., Aliu, E., Anderhub, H., et al. 2006, Science, 312, 1771.
  * Aliu et al. (2014) Aliu, E., Archambault, S., Aune, T., et al. 2014, ApJ, 780, 168.
  * Anderhub et al. (2009) Anderhub, H., Antonelli, L. A., Antoranz, P., et al. 2009, ApJ, 706, L27.
  * Aragona et al. (2009) Aragona, C., McSwain, M. V., Grundstrom, E. D., et al. 2009, ApJ, 698, 514.
  * Barkov & Bosch-Ramon (2018) Barkov, M. V., & Bosch-Ramon, V. 2018, MNRAS, 479, 1320.
  * Barkov & Bosch-Ramon (2021) Barkov, M. V., & Bosch-Ramon, V. 2021, Universe, 7, 277.
  * Bongiorno et al. (2011) Bongiorno, S. D., Falcone, A. D., Stroh, M., et al. 2011, ApJ, 737, L11.
  * Bosch-Ramon, Barkov, & Perucho (2015) Bosch-Ramon, V., Barkov, M. V., & Perucho, M. 2015, A&A, 577, A89.
  * Bosch-Ramon, Khangulyan, & Aharonian (2008) Bosch-Ramon, V., Khangulyan, D., & Aharonian, F. A. 2008, A&A, 482, 397.
  * Bosch-Ramon (2021) Bosch-Ramon, V. 2021, A&A, 645, A86.
  * Bosch-Ramon & Paredes (2004) Bosch-Ramon, V., & Paredes, J. M. 2004, A&A, 425, 1069.
  * Bosch-Ramon (2013) Bosch-Ramon, V. 2013, A&A, 560, A32.
  * Bosch-Ramon et al. (2012) Bosch-Ramon, V., Barkov, M. V., Khangulyan, D., & Perucho, M. 2012, A&A, 544, A59.
  * Bosch-Ramon et al. (2017) Bosch-Ramon, V., Barkov, M. V., Mignone, A., & Bordas, P. 2017, MNRAS, 471, L150.
  * Camilo et al. (2009) Camilo, F., Ray, P. S., Ransom, S. M., et al. 2009, ApJ, 705, 1.
  * Carciofi & Bjorkman (2006) Carciofi, A. C., & Bjorkman, J. E. 2006, ApJ, 639, 1081.
  * Carciofi & Bjorkman (2008) Carciofi, A. C., & Bjorkman, J. E. 2008, ApJ, 684, 1374.
  * Casares et al. (2005) Casares, J., Ribas, I., Paredes, J. M., Martí, J., & Allende Prieto, C. 2005, MNRAS, 360, 1105.
  * Casares et al. (2012) Casares, J., Ribó, M., Ribas, I., et al. 2012, MNRAS, 421, 1103.
  * Cañellas et al. (2012) Cañellas, A., Joshi, B. C., Paredes, J. M., et al. 2012, A&A, 543, A122.
  * Cerutti, Dubus, & Henri (2008) Cerutti, B., Dubus, G., & Henri, G. 2008, A&A, 488, 37.
  * Chen et al. (2021a) Chen, A. M., Guo, Y. D., Yu, Y. W., & Takata, J. 2021a, A&A, 652, A39.
  * Chen et al. (2021b) Chen, A. M., Ng, C. W., Takata, J., & Yu, Y. W. 2021b, RAA, 8, 189.
  * Chen et al. (2019) Chen, A. M., Takata, J., Yi, S. X., Yu, Y. W., & Cheng, K. S. 2019, A&A, 627, A87.
  * Chernyakova et al. (2006) Chernyakova, M., Neronov, A., Lutovinov, A., Rodriguez, J., & Johnston, S. 2006, MNRAS, 367, 1201.
  * Chernyakova et al. (2009) Chernyakova, M., Neronov, A., Aharonian, F., Uchiyama, Y., & Takahashi, T. 2009, MNRAS, 397, 2123.
  * Chernyakova et al. (2012) Chernyakova, M., Neronov, A., Molkov, S., et al. 2012, ApJ, 747, L29.
  * Chernyakova et al. (2014) Chernyakova, M., Abdo, A. A., Neronov, A., et al. 2014, MNRAS, 439, 432.
  * Chernyakova et al. (2015) Chernyakova, M., Neronov, A., van Soelen, B., et al. 2015, MNRAS, 454, 1358.
  * Chernyakova et al. (2017) Chernyakova, M., Babyk, I., Malyshev, D., et al. 2017, MNRAS, 470, 1718.
  * Chernyakova & Malyshev (2020) Chernyakova, M., & Malyshev, D. 2020, Multifrequency Behaviour of High Energy Cosmic Sources - XIII. 3-8 June 2019. Palermo, 45.
  * Chernyakova et al. (2020a) Chernyakova, M., Malyshev, D., Mc Keague, S., et al. 2020a, MNRAS, 497, 648.
  * Chernyakova et al. (2020b) Chernyakova, M., Malyshev, D., Blay, P., van Soelen, B., & Tsygankov, S. 2020b, MNRAS, 495, 365.
  * Chernyakova et al. (2021) Chernyakova, M., Malyshev, D., van Soelen, B., et al. 2021, Universe, 7, 242.
  * Coe et al. (2019) Coe, M. J., Okazaki, A. T., Steele, I. A., et al. 2019, MNRAS, 485, 1864.
  * Corbet et al. (2019) Corbet, R. H. D., Chomiuk, L., Coe, M. J., et al. 2019, ApJ, 884, 93.
  * Corbet et al. (2016) Corbet, R. H. D., Chomiuk, L., Coe, M. J., et al. 2016, ApJ, 829, 105.
  * de la Cita et al. (2017) de la Cita, V. M., Bosch-Ramon, V., Paredes-Fortuny, X., Khangulyan, D., & Perucho, M. 2017, A&A, 598, A13.
  * Dubus (2006a) Dubus, G. 2006a, A&A, 456, 801.
  * Dubus (2006b) Dubus, G. 2006b, A&A, 451, 9.
  * Dubus, Cerutti, & Henri (2010) Dubus, G., Cerutti, B., & Henri, G. 2010, A&A, 516, A18.
  * Dubus & Cerutti (2013) Dubus, G., & Cerutti, B. 2013, A&A, 557, A127.
  * Dubus (2013) Dubus, G. 2013, A&A Rev., 21, 64.
  * Dubus (2015) Dubus, G. 2015, Comptes Rendus Physique, 16, 661.
  * Dubus et al. (2017) Dubus, G., Guillard, N., Petrucci, P.-O., & Martin, P. 2017, A&A, 608, A59.
  * Eger et al. (2016) Eger, P., Laffon, H., Bordas, P., et al. 2016, MNRAS, 457, 1753.
  * Fermi LAT Collaboration et al. (2012) Fermi LAT Collaboration, Ackermann, M., Ajello, M., Ballet, et al. 2012, Science, 335, 189.
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., Prusti, et al. 2018, A&A, 616, A1.
  * Gregory & Taylor (1978) Gregory, P. C., & Taylor, A. R. 1978, Nature, 272, 704.
  * Gregory (2002) Gregory, P. C. 2002, ApJ, 575, 427.
  * Grundstrom et al. (2007) Grundstrom, E. D., Caballero-Nieves, S. M., Gies, D. R., et al. 2007, ApJ, 656, 437.
  * H. E. S. S. Collaboration et al. (2015a) H. E. S. S. Collaboration, Abramowski, A., Aharonian, F., Ait Benkhali, et al. 2015a, A&A, 577, A131.
  * H. E. S. S. Collaboration et al. (2015b) H. E. S. S. Collaboration, Abramowski, A., Acero, F., Aharonian, et al. 2015b, MNRAS, 446, 1163.
  * H. E. S. S. Collaboration et al. (2018) H. E. S. S. Collaboration, Abdalla, H., Abramowski, A., Aharonian, et al. 2018, A&A, 610, L17.
  * H. E. S. S. Collaboration et al. (2020) H. E. S. S. Collaboration, Abdalla, H., Adam, R., Aharonian, et al. 2020, A&A, 633, A102.
  * Harrison et al. (2000) Harrison, F. A., Ray, P. S., Leahy, D. A., Waltman, E. B., & Pooley, G. G. 2000, ApJ, 528, 454.
  * Hermsen et al. (1977) Hermsen, W., Swanenburg, B. N., Bignami, G. F., et al. 1977, Nature, 269, 494.
  * Hinton et al. (2009) Hinton, J. A., Skilton, J. L., Funk, S., et al. 2009, ApJ, 690, L101.
  * Ho et al. (2017) Ho, W. C. G., Ng, C.-Y., Lyne, A. G., et al. 2017, MNRAS, 464, 1211.
  * Huber et al. (2021a) Huber, D., Kissmann, R., Reimer, A., & Reimer, O. 2021a, A&A, 646, A91.
  * Huber, Kissmann, & Reimer (2021b) Huber, D., Kissmann, R., & Reimer, O. 2021b, A&A, 649, A71.
  * Hutchings & Crampton (1981) Hutchings, J. B., & Crampton, D. 1981, PASP, 93, 486.
  * Ignace & Brimeyer (2006) Ignace, R., & Brimeyer, A. 2006, MNRAS, 371, 343.
  * Jaron et al. (2018) Jaron, F., Massi, M., Kiehlmann, S., & Hovatta, T. 2018, MNRAS, 478, 440.
  * Jaron (2021) Jaron, F. 2021, Universe, 7, 245.
  * Johnston et al. (1994) Johnston, S., Manchester, R. N., Lyne, A. G., Nicastro, L., & Spyromilio, J. 1994, MNRAS, 268, 430.
  * Johnston et al. (1992) Johnston, S., Manchester, R. N., Lyne, A. G., et al. 1992, ApJ, 387, L37.
  * Kennel & Coroniti (1984a) Kennel, C. F., & Coroniti, F. V. 1984a, ApJ, 283, 694.
  * Kennel & Coroniti (1984b) Kennel, C. F., & Coroniti, F. V. 1984b, ApJ, 283, 710.
  * Kerschhaggl (2011) Kerschhaggl, M. 2011, A&A, 525, A80.
  * Khangulyan et al. (2007) Khangulyan, D., Hnatic, S., Aharonian, F., & Bogovalov, S. 2007, MNRAS, 380, 320.
  * Khangulyan, Aharonian, & Kelner (2014) Khangulyan, D., Aharonian, F. A., & Kelner, S. R. 2014, ApJ, 783, 100.
  * Khangulyan et al. (2012) Khangulyan, D., Aharonian, F. A., Bogovalov, S. V., & Ribó, M. 2012, ApJ, 752, L17.
  * Kirk, Ball, & Skjæraasen (1999) Kirk, J. G., Ball, L., & Skjæraasen, O. 1999, Astroparticle Physics, 10, 31.
  * Kong et al. (2011) Kong, S. W., Yu, Y. W., Huang, Y. F., & Cheng, K. S. 2011, MNRAS, 416, 1067.
  * Kong, Cheng, & Huang (2012) Kong, S. W., Cheng, K. S., & Huang, Y. F. 2012, ApJ, 753, 127.
  * Kravtsov et al. (2020) Kravtsov, V., Berdyugin, A. V., Piirola, V., et al. 2020, A&A, 643, A170.
  * Lamberts (2016) Lamberts, A. 2016, Bright Emissaries: Be Stars as Messengers of Star-Disk Physics, 506, 231.
  * Leahy (2004) Leahy, D. A. 2004, A&A, 413, 1019.
  * Li et al. (2011) Li, J., Torres, D. F., Zhang, S., et al. 2011, ApJ, 733, 89.
  * Li et al. (2018) Li, K. L., Takata, J., Ng, C. W., et al. 2018, ApJ, 857, 123.
  * Li et al. (2017) Li, K. L., Kong, A. K. H., Tam, P. H. T., et al. 2017, ApJ, 843, 85.
  * Lyne et al. (2015) Lyne, A. G., Stappers, B. W., Keith, M. J., et al. 2015, MNRAS, 451, 581.
  * Malyshev et al. (2019) Malyshev, D., Chernyakova, M., Santangelo, A., & Pühlhofer, G. 2019, Astronomische Nachrichten, 340, 465.
  * Maraschi & Treves (1981) Maraschi, L., & Treves, A. 1981, MNRAS, 194, 1P.
  * Marcote (2017) Marcote, B. 2017, XII Multifrequency Behaviour of High Energy Cosmic Sources Workshop (MULTIF2017), 45.
  * Massi, Ros, & Zimmermann (2012) Massi, M., Ros, E., & Zimmermann, L. 2012, A&A, 540, A142.
  * Massi & Torricelli-Ciamponi (2016) Massi, M., & Torricelli-Ciamponi, G. 2016, A&A, 585, A123.
  * Massi & Torricelli-Ciamponi (2014) Massi, M., & Torricelli-Ciamponi, G. 2014, A&A, 564, A23.
  * McSwain et al. (2011) McSwain, M. V., Ray, P. S., Ransom, S. M., et al. 2011, ApJ, 738, 105.
  * Miller-Jones et al. (2018) Miller-Jones, J. C. A., Deller, A. T., Shannon, R. M., et al. 2018, MNRAS, 479, 4849.
  * Moderski et al. (2005) Moderski, R., Sikora, M., Coppi, P. S., & Aharonian, F. 2005, MNRAS, 363, 954.
  * Molina & Bosch-Ramon (2020) Molina, E., & Bosch-Ramon, V. 2020, A&A, 641, A84.
  * Moritani et al. (2018) Moritani, Y., Kawano, T., Chimasu, S., et al. 2018, PASJ, 70, 61.
  * Moritani & Kawachi (2021) Moritani, Y., Kawachi, A. 2021, Universe, 7, 320.
  * Motch et al. (1997) Motch, C., Haberl, F., Dennerl, K., Pakull, M., & Janot-Pacheco, E. 1997, A&A, 323, 853.
  * Negueruela et al. (2011) Negueruela, I., Ribó, M., Herrero, A., et al. 2011, ApJ, 732, L11.
  * Ng et al. (2019) Ng, C.-Y., Ho, W. C. G., Gotthelf, E. V., et al. 2019, ApJ, 880, 147.
  * Pal et al. (2019) Pal, P. S., Tam, P. H. T., Cui, Y., et al. 2019, ApJ, 882, 25.
  * Papitto, Torres, & Rea (2012) Papitto, A., Torres, D. F., & Rea, N. 2012, ApJ, 756, 188.
  * Paredes-Fortuny et al. (2015) Paredes-Fortuny, X., Bosch-Ramon, V., Perucho, M., & Ribó, M. 2015, A&A, 574, A77.
  * Paredes & Figueras (1986) Paredes, J. M., & Figueras, F. 1986, A&A, 154, L30.
  * Paredes & Bordas (2019) Paredes, J. M., & Bordas, P. 2019, Frontier Research in Astrophysics - III. 28 May - 2 June 2018. Mondello Palermo, 44.
  * Petropoulou et al. (2018) Petropoulou, M., Vasilopoulos, G., Christie, I. M., Giannios, D., & Coe, M. J. 2018, MNRAS, 474, L22.
  * Porter & Rivinius (2003) Porter, J. M., & Rivinius, T. 2003, PASP, 115, 1153.
  * Rivinius, Carciofi, & Martayan (2013) Rivinius, T., Carciofi, A. C., & Martayan, C. 2013, A&A Rev., 21, 69.
  * Romoli et al. (2015) Romoli, C., Bordas, P., Mariaud, C., et al. 2015, 34th International Cosmic Ray Conference (ICRC2015), 34, 873.
  * Saha et al. (2016) Saha, L., Chitnis, V. R., Shukla, A., Rao, A. R., & Acharya, B. S. 2016, ApJ, 823, 134.
  * Shannon, Johnston, & Manchester (2014) Shannon, R. M., Johnston, S., & Manchester, R. N. 2014, MNRAS, 437, 3255.
  * Sidoli et al. (2006) Sidoli, L., Pellizzoni, A., Vercellone, S., et al. 2006, A&A, 459, 901.
  * Takata & Taam (2009) Takata, J., & Taam, R. E. 2009, ApJ, 702, 100.
  * Takata et al. (2012) Takata, J., Okazaki, A. T., Nagataki, S., et al. 2012, ApJ, 750, 70.
  * Takata et al. (2014) Takata, J., Leung, G. C. K., Tam, P. H. T., et al. 2014, ApJ, 790, 18.
  * Takata et al. (2017) Takata, J., Tam, P. H. T., Ng, C. W., et al. 2017, ApJ, 836, 241.
  * Tam et al. (2018) Tam, P. H. T., He, X.-B., Pal, P. S., & Cui, Y. 2018, ApJ, 862, 165.
  * Tavani & Arons (1997) Tavani, M., & Arons, J. 1997, ApJ, 477, 439.
  * Tokayer et al. (2021) Tokayer, Y. M., An, H., Halpern, J. P., et al. 2021, arXiv e-prints, arXiv:2110.01075.
  * Torres et al. (2012) Torres, D. F., Rea, N., Esposito, P., et al. 2012, ApJ, 744, 106.
  * van Soelen (2017) van Soelen, B. 2017, 5th Annual Conference on High Energy Astrophysics in Southern Africa, 31.
  * van Soelen & Meintjes (2011) van Soelen, B., & Meintjes, P. J. 2011, MNRAS, 412, 1721.
  * van Soelen et al. (2012) van Soelen, B., Meintjes, P. J., Odendaal, A., & Townsend, L. J. 2012, MNRAS, 426, 3135.
  * Wang et al. (2004) Wang, N., Johnston, S., & Manchester, R. N. 2004, MNRAS, 351, 599.
  * Weng et al. (2021) Weng, S.-S., Pan, Z., Qian, L., et al. 2021, The Astronomer’s Telegram, 14297, 1.
  * Xing, Wang, & Takata (2017) Xing, Y., Wang, Z., & Takata, J. 2017, ApJ, 851, 92.
  * Yi & Cheng (2017) Yi, S.-X., & Cheng, K. S. 2017, MNRAS, 471, 4228.
  * Zabalza, Paredes, & Bosch-Ramon (2011) Zabalza, V., Paredes, J. M., & Bosch-Ramon, V. 2011a, A&A, 527, A9.
  * Zabalza, Bosch-Ramon, & Paredes (2011) Zabalza, V., Bosch-Ramon, V., & Paredes, J. M. 2011b, ApJ, 743, 7.
  * Zdziarski, Neronov, & Chernyakova (2010) Zdziarski, A. A., Neronov, A., & Chernyakova, M. 2010, MNRAS, 403, 1873.
  * Zhang, Zhang, & Fang (2009) Zhang, J.-F., Zhang, L., & Fang, J. 2009, PASJ, 61, 949.