# Two models for the orbital modulation of $\gamma$-rays in Cyg X-3

Anton Dmytriiev Centre for Space Research, North-West University,
Potchefstroom, 2520, South Africa<EMAIL_ADDRESS>Andrzej A.
Zdziarski Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences,
Bartycka 18, PL-00-716 Warszawa, Poland<EMAIL_ADDRESS>Denys Malyshev
Institut für Astronomie und Astrophysik Tübingen, Universität Tübingen, Sand
1, D-72076 Tübingen, Germany  Valentí Bosch-Ramon Departament de Física
Quàntica i Astrofísica, Institut de Ciències del Cosmos, Universitat de
Barcelona, Martí i Franquès 1, 08028 Barcelona, Spain  Instituto Argentino de
Radioastronomía, CONICET, Berazategui Partido, Buenos Aires Province,
Argentina  Maria Chernyakova School of Physical Sciences and Centre for
Astrophysics & Relativity, Dublin City University, Glasnevin, D09 W6Y4,
Ireland  Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin
2, Ireland

###### Abstract

We model the currently available $\gamma$-ray data from the Fermi Large Area
Telescope on Cyg X-3. Thanks to its very strong $\gamma$-ray activity during
2018–2021, the data quality has significantly improved. We study the strong
orbital modulation of the $\gamma$-rays observed during at high $\gamma$-ray
fluxes. The modulation, as found earlier, is well modeled by anisotropic
Compton scattering of the donor blackbody emission by relativistic electrons
in a jet strongly misaligned with respect to the orbital axis. We confirm that
this model fits well both the average $\gamma$-ray modulation light curve and
the spectrum. However, we find that if the jet is aligned with the spin axis
of a rotating black hole, it would undergo geodetic precession with the period
of $\sim$50 years. However, its presence is ruled out by both the $\gamma$-ray
and radio data. Therefore, we consider an alternative model in which the
average jet direction jet is aligned, but it is bent to outside the orbit
owing to the thrust of the donor stellar wind, and thus precesses at the
orbital period. The $\gamma$-ray modulation appears then owing to the variable
Doppler boosting of synchrotron self-Compton jet emission. The model also fits
well the data. However, the fitted bending angle is much larger than the
theoretical one based on the binary and wind parameters as currently known.
Thus, both models disagree with important aspects of our current theoretical
understanding of the system. We discuss possible ways to find the correct
model.

revtex4-1.clsrevtex4-2.cls

##  1 Introduction

The high-mass binary Cyg X-3 is one of the first discovered (Giacconi et al.,
1967) X-ray binaries. Still, after many years of intense research, it remains
to be a puzzling (and unique) system. The nature of its compact object remains
uncertain. It can be either a neutron star or a black hole (BH), see Koljonen
& Maccarone (2017) and Antokhin et al. (2022), hereafter A22, for discussions.
Still, those authors favored the presence of a BH, which is also supported by
various aspects of the X-ray and radio emission (Hjalmarsdotter et al., 2008,
2009; Szostek & Zdziarski, 2008; Szostek et al., 2008; Koljonen et al., 2010).
Its donor is a Wolf-Rayet (WR) star (van Kerkwijk et al., 1992, 1996; Koljonen
& Maccarone, 2017), which makes it the only known binary system in the Milky
Way consisting of a WR star and a compact object, and a candidate for a future
merger (Belczyński et al., 2013). The compact object accretes a fraction of
the stellar wind from the donor. A22 estimated the total mass of the system as
$18.8{M}_{\sun}$, and the compact object mass as $M_{\rm c}\approx
7.2{M}_{\sun}$, implying it is a BH. However, they also considered lower
masses possible if the wind is significantly clumped. Its distance has
recently been determined based on a radio parallax as $D\approx 9.7\pm 0.5$
kpc and based on Galactic proper motions and line-of-sight radial velocity
measurements as $9\pm 1$ kpc (Reid & Miller-Jones, 2023). We hereafter assume
$D=9$ kpc, $M_{*}=11.6{M}_{\sun}$, $M_{\rm c}=7.2{M}_{\sun}$.

Cyg X-3 is also the brightest and most variable radio source among X-ray
binaries, showing resolved relativistic jets (Mioduszewski et al., 2001;
Miller-Jones et al., 2004). Its jets also emit high-energy $\gamma$-rays, as
first discovered by Mori et al. (1997), and later confirmed by observations
with the Fermi Large Area Telescope (LAT; Fermi LAT Collaboration et al.
2009), and AGILE (Tavani et al., 2009). As found by Fermi LAT Collaboration et
al. (2009), the $\gamma$-rays emitted during flaring epochs are strongly
orbitally modulated. The modulated emission was explained by Dubus et al.
(2010) as anisotropic Compton scattering of blackbody photons from the donor
by relativistic electrons with a power-law distribution accelerated in the jet
(which we hereafter refer to as either the blackbody Compton model or Model
1). The peak of the emission of a jet aligned with the orbital axis would be
at the superior conjunction of the compact object. However, Dubus et al.
(2010) found that the peak was clearly before the superior conjunction,
indicating that the jet is inclined with respect to the orbital axis. With the
number of high-energy $\gamma$-rays increasing with time, the statistical
accuracy of the orbital modulation increased, and this model was fitted to the
updated data by Zdziarski et al. (2018), hereafter Z18.

Figure 1: The LAT high-energy $\gamma$-ray light curve of Cyg X-3 from the
beginning of the Fermi observations until MJD 60200. The red and blue symbols
represent the detections within a day and upper limits, respectively. Our
criterion defining the $\gamma$-ray bright state for 1-d integration is shown
by the horizontal dotted line. Z18 analyzed the data until MJD 57982, which is
marked by the vertical dotted line. We see that the later observations have
multiplied the number of detected photons by a substantial factor.

The observations analyzed by Z18 were until MJD 57982. As we have found out by
analyzing the LAT data, Cyg X-3 has shown very strong activity after that
date, see Figure 1. The new data allow us to increase the number of detected
$\gamma$-ray photons by a large factor, and, consequently, to strongly reduce
the errors on the light curve averaged and folded over the orbital period.

In our modeling, we consider an important consequence of the jet–binary axis
misalignment. Jets powered by the rotation of BHs (Blandford & Znajek, 1977)
are launched along the BH spin axis (McKinney et al., 2013). Jets launched via
poloidal magnetic field of accretion disks (Blandford & Payne, 1982) will
initially be directed along the inner disk rotation axis, which could, in
turn, be also aligned with the BH spin (Bardeen & Petterson, 1975). In those
cases, the observed jet–orbit misalignment would correspond to the BH
spin–orbit one. However, spin precession due to the binary angular
momentum–compact object spin coupling has to take place in general relativity,
and the period of such geodetic precession for Cyg X-3 is as short as $\sim$50
years. Such precession is clearly not seen in the Fermi data, which span
$\approx$14 years, i.e., more than a quarter of that period. It is also not
seen in the $\sim$30 years of the radio imaging data, in which the jet
position angle remains constant along the north-to-south direction. As shown
in this work, these observations present an argument against the misaligned
jet model.

Therefore, we also consider an alternative model of the orbital modulation in
which it is due to the Doppler boosting of the jet synchrotron self-Compton
(SSC) emission. This could be achieved if the jet is bent to outside by the
thrust of the stellar wind (Yoon & Heinz, 2015; Yoon et al., 2016; Bosch-Ramon
& Barkov, 2016; Molina & Bosch-Ramon, 2018; Molina et al., 2019; Barkov &
Bosch-Ramon, 2022), with the wind being very strong in this binary. In
addition, the Coriolis force due to the binary rotation bends the jet opposite
to the rotation direction. We hereafter refer to it as either the orbital
precession model or Model 2. We have found that the resulting modulated
Doppler boosting can fit the data at an accuracy similar to the blackbody
Compton model. In this model, the average jet direction is fully aligned with
the orbital axis, and thus no geodetic precession is expected. An important
prediction of this model is an orbital modulation of the associated
synchrotron emission, with the orbital light curve similar to that of the
$\gamma$-rays. However, we find that the fitted bending angle is much larger
than that expected theoretically based on the likely wind and jet parameters.
This, in turn, represents a strong argument against this model.

We present our detailed results for both models. We then discuss possible
solutions of their problems for each of them. According to our present
knowledge, we cannot determine which model is more likely.

Figure 2: The Lomb-Scargle periodogram in the $\gamma$-ray bright state in the
0.1–100 GeV range, calculated accounting for the secular orbital period
increase and taking into account the measurement uncertainties. The only
significant peak is equal within the fitted uncertainties to the orbital
period of Cyg X-3.

##  2 The data analysis

We analyze the data from the direction of Cyg X-3 for MJD 57982–60201, see
Figure 1. We used the Fermi Science Tools v. 2.2.0. Details of the analysis
generally follow those given in Z18. For the spectral analysis, we used the
energy range of 0.05–500 GeV and for the timing analysis, that of 0.1–100 GeV.
In this paper, we analyze only $\gamma$-ray bright states, which we define in
the same way as the ‘flaring state’ in Z18. The minimum level of this state is
shown by the horizontal line in Figure 1. These fluxes correspond to
$F(0.1$–100 GeV$)>5\times 10^{-10}$ erg cm-2 s-1 and the test statistics
(Mattox et al., 1996) of ${\rm TS}>16$, i.e., the significance $>4\sigma$. We
obtain the average spectrum for that state with positive detections in the
0.05–20 GeV range, and upper limits at higher energies.

We used the quadratic ephemeris given by model 4 of Antokhin & Cherepashchuk
(2019) and searched for periodicity in the LAT light curve taking into account
their rate of the period increase. We used the Lomb-Scargle method. Our
results are shown in Figure 2. We have found the period of
$P_{0}=0.1996847(2)$ d, which is in full agreement with that of Antokhin &
Cherepashchuk (2019), $P_{0}=0.199684622(15)$ d. We also tested for the effect
of using their quadratic+sinusoidal ephemeris, but we found its effect to be
very minor. Then, we calculated the folded and averaged light curve in 10
phase bins.

##  3 Theoretical description

Figure 3: The geometry corresponding to emission in a compact region of the
jet (the counterjet is not shown) for the assumed counter-clockwise rotation.
The axes $x$ and $y$ are in the binary plane, and the $+z$ direction is
parallel to the orbital axis. The $-x$ direction gives the projection of the
direction toward the observer onto the binary plane. The observer is at an
angle, $i$, with respect to the orbital axis, $\varphi$ is the orbital phase
(the superior conjunction is at $\varphi=0$, and the shown configuration is
close to the inferior conjunction), $\theta_{\rm j}$ and $\varphi_{\rm j}$ are
the polar and azimuthal angles of the jet with respect to the orbital axis,
$a$ is the binary separation, $h$ is the distance of the $\gamma$-ray source
from the center of the compact object, ${\mathbfit e}_{\rm orb}$ is a unit
vector perpendicular to the orbital plane, ${\mathbfit e}_{\rm obs}$,
${\mathbfit e}_{\rm c}$ and ${\mathbfit e}_{*}$ point from the donor toward
the observer, the center of the compact object, and the $\gamma$-ray source,
respectively, and ${\mathbfit e}_{\rm j}$ points from the center of the
compact object toward the $\gamma$-ray source.

###  3.1 The geometry

Our assumed geometry is shown in Figure 3. The binary rotates with the period
$P$, and the azimuth of the rotation with respect to the superior conjunction
is $\varphi$. We assume no eccentricity and the rotation to be counter-
clockwise, as found by Veledina et al. (2024). The $\gamma$-ray emitting
region in the jet is at the distance $h$ from the compact object, which itself
is at the binary separation, $a$, from the center of the donor. The jet is
inclined with respect to the orbital axis as given by the polar and azimuthal
angles, $\theta_{\rm j}$ and $\varphi_{\rm j}$, respectively, and the distance
of the emitting region from the center of the donor is $r$. The jet polar
angle, $\varphi_{\rm j}$, is either fixed in the case of the blackbody Compton
model or dependent on the orbital phase in the orbital precession model.

The unit vectors pointing towards the observer, from the stellar center to the
compact object, along the jet, and along $z$ (parallel to the orbital axis),
are

$\displaystyle{\mathbfit e}_{\rm obs}=(-\sin i,0,\cos i),\quad{\mathbfit
e}_{\rm c}=(\cos\varphi,\sin\varphi,0),$ (1) $\displaystyle{\mathbfit e}_{\rm
j}=(\sin\theta_{\rm j}\cos\varphi_{\rm j},\sin\theta_{\rm j}\sin\varphi_{\rm
j},\cos\theta_{\rm j}),\quad{\mathbfit e}_{\rm orb}=(0,0,1),$ (2)

respectively. We make a simplifying assumption that the modulated emission at
a given frequency originates at a single distance, $h$, from the BH. The
vector connecting the donor center with the emission point is $a{\mathbfit
e}_{\rm c}+h{\mathbfit e}_{\rm j}$. The length of this vector gives the
distance of that point from the donor, and its square is given by

$r^{2}=h^{2}+a^{2}+2ha\sin\theta_{\rm j}\cos(\varphi_{\rm j}-\varphi).$ (3)

An important issue concerns the determination of the phase of the superior
conjunction. All of the determinations of the ephemeris so far have been based
on fitting average profiles of the X-ray emission in various bands (see
Bhargava et al. 2017; Antokhin & Cherepashchuk 2019 for the most recent work).
The X-ray modulation appears due to the bound-free absorption and scattering
by the stellar wind emitted by the donor. Since those profiles are not
symmetric, a template designed by van der Klis & Bonnet-Bidaud (1989) was used
for the fits. However, the phases corresponding to the observed minimum fluxes
using that template are not null (see, e.g., the profiles observed in various
X-ray bands and the hard and soft X-ray states of the source in Zdziarski et
al. 2012b). This issue was studied in detail by A22. They used the X-ray data
from the All Sky Monitor (Levine et al., 1996) onboard Rossi X-ray Timing
Explorer, and the infrared data acquired by them. They estimated the superior
conjunction based on their three-component wind model fitted to those data to
be $\phi_{0}/2\pi=-0.066\pm 0.006$, where $\phi$ is the orbital phase defined
by the ephemeris. Since the highest $\gamma$-ray flux in our folded/averaged
light curve is in the phase bin with $\phi/2\pi=-0.15\pm 0.05$, the correction
to the phase of the superior conjunction is significant. Thus, the true
orbital phase of the binary is given by

$\varphi/2\pi=\phi/2\pi+0.066\pm 0.006.$ (4)

When the jet is inclined with respect to the orbital axis, its viewing angle
is different from that of the orbital axis itself. It equals to

$i_{\rm j}=\arccos(\cos i\cos\theta_{\rm j}-\cos\varphi_{\rm j}\sin
i\sin\theta_{\rm j}).$ (5)

Furthermore, the position angle on the sky of the binary differs from that of
the jet. The difference between the jet and the binary position angles,
$\Delta\lambda$, is given by

$\cos\Delta\lambda=\frac{\cos i\cos\varphi_{\rm j}\sin\theta_{\rm j}+\sin
i\cos\theta_{\rm j}}{\sqrt{1-\left(\cos i\cos\theta_{\rm j}-\sin
i\cos\varphi_{\rm j}\sin\theta_{\rm j}\right)^{2}}}.$ (6)

The value of the binary separation follows via the Kepler law from the assumed
total mass of $18.8{M}_{\sun}$, yielding $a\approx 2.66\times 10^{11}$ cm. We
either assume the binary inclination of $i=30\arcdeg$ (A22) or treat it as a
free parameter.

###  3.2 The electrons

The electrons in the emitting regions are assumed to be isotropic in the
comoving frame, and move with bulk the velocity $\beta_{\rm j}c$. Also, they
are assumed to have a power law distribution, $N(\gamma)=K\gamma^{-p}$ from
$\gamma_{\rm min}$ to $\gamma_{\rm max}$. Hereafter, $N(\gamma)$ refers to all
of the power-law electrons in the source. However, we normalize the electron
distribution to the kinetic energy content (i.e. excluding the rest energy)
between $\gamma=1$ to $\gamma_{\rm max}$, $E_{\rm e}$, instead of $K$.

We also take into account the electrons cooled by the blackbody Compton,
synchrotron, SSC and adiabatic losses to below $\gamma_{\rm min}$. This is
important since these particles contribute significantly to the emission at
low energies. The electron distribution below $\gamma_{\rm min}$ is determined
by a kinetic equation, where the radiative and adiabatic losses compete.
Solving the kinetic equation with no injection of fresh electrons below
$\gamma_{\rm min}$ yields approximately a broken power-law distribution, with
$N_{\rm e}\propto\gamma^{-2}$ in the range $\gamma_{\rm
br}\leq\gamma<\gamma_{\rm min}$ (where the radiative losses dominate), and
$N_{\rm e}\propto\gamma^{-1}$ in the range $1\lesssim\gamma<\gamma_{\rm br}$
(where the adiabatic losses dominate). The former follows from both the
Compton losses in the Thomson regime and the synchrotron ones having
$\dot{\gamma}\propto-\gamma^{2}$. The position of the cooling break
$\gamma_{\rm br}$ is determined self-consistently in each simulation by
equating the timescales of total radiative losses and adiabatic losses. We
verify that the Compton scattering around $\gamma_{\rm br}$ is in the Thomson
regime.

###  3.3 Radiation

We assume a one-zone $\gamma$-ray emitting region in our model. To validate
this assumption, we examined the importance of propagation effects discussed
in Khangulyan et al. (2018), which, if significant, might lead to a departure
from the one-zone model. We find that these propagation effects have only a
minor impact in our specific case. This is due to the fact that the jet
velocity $\beta_{\rm j}$ in all our fits remains below the value identified in
Khangulyan et al. (2018), above which such effects become non-negligible.

The treatment of the anisotropic Compton scattering off stellar blackbody
photons follows that of Zdziarski et al. (2012c), which employs the full
Klein-Nishina (KN) cross section integrated over a broken power law electron
distribution. The blackbody photons are calculated for a donor with radius
$R_{*}$ and temperature $T_{*}$, but we approximate that emission as coming
from a point source. For that process, the total number of electrons
determines the observed flux. The emission region has to have a low vertical
extent since the orbital profile of the blackbody Compton emission depends
sensitively on the height.

For an assumed value of the magnetic field strength, $B$, we compute
synchrotron and SSC emission assuming that the jet emitting region has a
cylindrical geometry and is filled with a tangled and uniform magnetic field.
The emission is therefore isotropic in the comoving frame. The radius of the
cylinder is defined by the opening angle of the jet $\alpha_{\rm j}$ as
$R=\alpha_{\rm j}h$, and for the height of the cylinder we arbitrarily choose
$\Delta h=(1/3)h$. We note that the observed synchrotron spectrum for given
$B$ and $K$ does not depend on neither $\Delta h$ nor $h$, with its luminosity
being $L_{\rm S}\propto B^{2}K$, apart from the synchrotron self-absorbed
domain at low frequencies. Then, the shape of the SSC spectrum is also
independent of them. However, the ratio of the two luminosities (the Compton
dominance), is

$\frac{L_{\rm SSC}}{L_{\rm S}}\propto\frac{L_{\rm S}}{\alpha_{\rm j}h\Delta
hB^{2}}\propto\frac{K}{\alpha_{\rm j}h\Delta h}.$ (7)

In our Model 2, with precessing jet, $L_{\rm SSC}$ is fixed by the observed
$\gamma$-ray flux, and hence the predicted $L_{\rm S}$ increases when
increasing $\alpha_{\rm j}h\Delta h$.

For the single electron synchrotron emissivity, we employ an approximation
from eq. (3.40) in Boettcher et al. (2012). We include the synchrotron self-
absorption effect by using eq. (3.42) in Boettcher et al. (2012). To obtain
the synchrotron radiation energy density, we use the solution of the radiation
transfer equation for cylindrical geometry (e.g., eq. 14 in Chiaberge &
Ghisellini, 1999), while to calculate the observed synchrotron flux in the
absorbed regime, we follow eqs. (18–19) of Zdziarski et al. (2012a). The SSC
computation takes into account the KN cross section and follows eqs. (8–12) in
Katarzyński et al. (2001) (though excluding the scaling factor of 3/4 for
spherical geometry present in their eq. 11). The emission is transformed to
the observer’s frame, including contributions from both the jet and
counterjet.

We also check the importance of internal $\gamma$–$\gamma$absorption of high-
energy $\gamma$-rays on the soft synchrotron photon field. We find that this
effect is negligible within our studied energy range, and becomes pronounced
only above 10–20 GeV. At those energies, however, the the source is no longer
detected.

###  3.4 Additional observational constraints

Apart from the LAT $\gamma$-ray data, we impose some observational constraints
in other spectral regions. First, the X-ray $EF_{E}$ flux at 100 keV measured
in the soft states of Cyg X-3 is $\approx$0.2 keV cm-2 s-1 (Szostek et al.,
2008). The observed $\sim$100 keV flux exhibits strong orbital modulation
(Zdziarski et al., 2012b) similar to that at lower energies, with the minimum
around the superior conjunction, which could be due to wind absorption. This
implies that most of the X-ray emission originates in the accretion flow
rather than in the jet, whose emission shows a very different modulation
pattern. On the other hand, the duration of the flaring state is by a factor
of $\sim$7 shorter than the one of the soft state. Thus, we require that jet
contribution at 100 keV to be $<$0.1 keV cm-2 s${}^{-1}\approx 1.6\times
10^{-10}$ erg cm-2 s-1.

Furthermore, the model synchrotron flux in the flaring state cannot exceed
available measurements. We use the observational constraint from the
Submillimeter Array (SMA) observations during flaring periods in 2018, with
the average flux of $\sim$300 mJy at 225 GHz (McCollough, 2023). That emission
undergoes free-free absorption on electrons of the stellar wind. This process
can influence both the observed level of 225 GHz emission and the profile of
its orbital modulation. We estimate the distance away from the star at which
the optical depth of the free-free absorption of 225 GHz emission drops to
unity. Using typical values for the mass loss rate of the donor, velocity and
temperature of the stellar wind in Cyg X-3 (e.g., A22), and assuming a fully
ionized helium plasma, we find this critical distance to be $\sim 20a$. This
provides an approximate constraint on the location of the region responsible
for producing the 225 GHz emission.

###  3.5 The models

We consider two models.

Model 1. The jet is assumed to have a fixed orientation. The dominant
radiative process is the anisotropic Compton scattering of stellar blackbody
photons in a compact region along the jet. The orbital modulation is entirely
due to this process since the Doppler boosting does not depend on the orbital
phase. As found before (Dubus et al. 2010; Z18), this model fits well the
$\gamma$-ray folded light curve, while requiring the jet to be significantly
misaligned with respect to the orbital axis. Since the jet also contains
magnetic field, the SSC contributes to the folded light curve at a constant
level, which effect was not calculated before (except for a qualitative
consideration in Zdziarski et al. 2012c). Here, we calculate the maximum
allowed field strength allowed by the data. The main parameters of this model
are the jet angles and the height of the emission region along the jet. The
emission of both the jet and the counterjet are taken into account, and we
have

$\varphi_{\rm j}={\rm const},\quad{\mathbfit e}_{\rm cj}=-{\mathbfit e}_{\rm
j},\quad\varphi_{\rm cj}=\pi+\varphi_{\rm j},\quad\theta_{\rm
cj}=\pi-\theta_{\rm j}.$ (8)

Model 2. The jet average direction is aligned with the orbital axis but the
thrust of the stellar wind bends the jet outward by an angle $\theta_{\rm j}$
along the line connecting the jet with the center of the donor (Yoon & Heinz,
2015; Yoon et al., 2016; Bosch-Ramon & Barkov, 2016; Barkov & Bosch-Ramon,
2022). This effect alone would result in a jet precessing at the orbital
period with the azimuth given by $\varphi_{\rm j}=\varphi$. However, the
Coriolis force arising from the orbital motion will induce a lateral bending
of the jet opposite to the direction of the rotation, leading to formation of
a helix. Here, we are considering emission of a compact region of the jet,
which can be then offset by an arbitrary angle, $\varphi_{\rm o}$. We note
that the Coriolis force acts in the same direction on the jet and counterjet,
which is evident upon considering the relevant vector product. This implies
that both jet and counterjet experience lateral bending opposite to the
direction of rotation. Then, we have

$\displaystyle\varphi_{\rm j}=\varphi-\varphi_{\rm o},$ (9)
$\displaystyle{\mathbfit e}_{\rm cj}=(\sin\theta_{\rm j}\cos\varphi_{\rm
j},\sin\theta_{\rm j}\sin\varphi_{\rm j},-\cos\theta_{\rm j}),$ (10)
$\displaystyle\varphi_{\rm cj}=\varphi_{\rm j},\quad\theta_{\rm
cj}=\pi-\theta_{\rm j}.$ (11)

The variable $\varphi_{\rm j,cj}$ results then in a variability of the Doppler
boosting, which is then responsible for the orbital modulation. We note that
our treatment of the geometry is approximate, since we assume a straight jet
connecting the origin with the emission region, while the jet actually has a
helical shape. The exact details of this helical shape are influenced by the
jet acceleration profile, i.e., how the jet acceleration rate (and hence the
jet velocity) evolves along the length of the jet. Since we do not have
precise information on this acceleration profile, we have to neglect this
complication.

A recent study by Yang et al. (2023) suggests the presence of a non-negligible
transversal component in the jet of Cyg X-3, albeit observed in the hard
state. Nevertheless, it is possible that this transversal component may arise
from the wind-induced jet precession, and hence the apparently broadened jet
in those observations may be attributed to its helical nature, supporting our
assumptions about the jet morphology in the Model 2.

The jet will still be irradiated by the stellar blackbody, which adds a
perturbation to the modulation due to the variable Doppler boosting of the
precessing jet. We calculate the maximum flux, $\propto R_{*}^{2}T_{*}^{4}$,
of the blackbody allowed by the model for a given $h$, constraining in turn
the radius and temperature of the donor. The maximum depends on the height of
the emission region, $h$.

We assess the significance of the attenuation of $\gamma$-rays en route to the
observer as they traverse the photon field of the donor star and suffer
$\gamma$–$\gamma$absorption. Our analysis shows that this effect is entirely
negligible for the jet, with substantial opacity emerging only beyond
$\sim$100 GeV. For the counterjet, the onset of significant attenuation is,
however, found to arise already at 5–6 GeV for Model 1, and 10–15 GeV for
Model 2. Still, considering its subdominant contribution to the total emission
and the increasingly large uncertainties in the data beyond these energies, we
conclude that, in our specific case, this effect can be neglected also for the
counterjet.

With either model, we conduct fitting of the phase-averaged spectrum and the
light curve. We first obtain the spectrum of the emission as a function of the
orbital phase for assumed preliminary values of $\theta_{\rm j}$, either
$\varphi_{\rm j}$ or $\varphi_{\rm o}$, $\beta_{\rm j}$, $i$, and $h$. With
that, we fit the observed phase-averaged photon spectrum, obtaining $E_{\rm
e}$, $p$, $\gamma_{\rm min}$ and $\gamma_{\rm max}$. With those, we fit in
turn the observed folded light curve, which yields updated values of the
parameters. Using those more accurate estimates, we fit again the phase-
averaged spectrum, iterating in this way until convergence. The fitting is
performed using the python LMFIT module (employing the default Levenberg-
Marquardt minimization method).

##  4 Results for Model 1: Blackbody Compton

###  4.1 Fit results

Table 1: Fit results for the phase-averaged spectrum and the light curve for Model 1 without and with magnetic field  $\alpha_{\rm j}(\arcdeg)$ |  $B_{\rm max}$(G) | $h/a$ | $i(\arcdeg)$ | $\theta_{\rm j}(\arcdeg)$ | $\varphi_{\rm j}(\arcdeg)$ | $\beta_{\rm j}$ | $\chi_{\nu}^{2}$ | $p$ | $\gamma_{\rm min}\,(10^{3})$ |  $E_{\rm e}\,(10^{37}$erg) | $\chi_{\nu}^{2}$ | $\gamma_{\rm br}$ | $i_{\rm j}(\arcdeg)$ | $\Delta\lambda(\arcdeg)$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
– | 0 | $1.96\pm 0.09$ | $28\pm 6$ | $32\pm 7$ | $188\pm 2$ | $0.59\pm 0.05$ | 7.2/5 | $4.16\pm 0.05$ | $2.9\pm 0.2$ | $18\pm 1$ | 5.1/7 | 426 | $5_{-3}^{+12}$ | $228_{-35}^{+118}$
5 | 34 | $2.3\pm 0.6$ | $33\pm 7$ | $35\pm 8$ | $188\pm 3$ | $0.7\pm 0.2$ | 7.3/5 | $4.8\pm 0.1$ | $3.4\pm 0.4$ | $14\pm 2$ | 7.5/7 | 591 | $5_{-2}^{+12}$ | $261_{-67}^{+88}$

\begin{overpic}[height=138.76157pt]{fit_SED_CygX3_Stable_Jet_zeroB.jpg}
\put(0.08,0.285){\large(a)} \end{overpic}

\begin{overpic}[height=138.76157pt]{fit_LC_CygX3_Stable_Jet_zeroB.jpg}
\put(0.41,0.28){\large(b)} \end{overpic}

\begin{overpic}[height=138.76157pt]{fit_SED_CygX3_Stable_Jet_maxB_34G_5deg.jpg}
\put(0.08,0.285){\large(c)} \end{overpic}

\begin{overpic}[height=138.76157pt]{fit_LC_CygX3_Stable_Jet_maxB_34G_5deg.jpg}
\put(0.41,0.28){\large(d)} \end{overpic}

Figure 4: Fit results for Model 1 without and with magnetic field. Black error
bars show the LAT data, the purple SED data point (upper limit) shows the
constraint on the X-ray flux at 100 keV of $\approx 1.6\times 10^{-10}$ erg
cm-2 s-1, the dark blue solid curves show the best-fit models, and the dashed
histograms show the model predictions for the binned fluxes. Cyan dashed and
the grey dash-dotted lines display the contributions from the jet and
counterjet, respectively. (a) The phase-averaged spectrum for $B=0$. (b) The
0.1–100 GeV light curve folded and averaged over the orbital phase with
respect to the superior conjunction, $\varphi$. (c) Similar to (a), but for
the upper limit value of the magnetic field $B_{\rm max}=34$ G. (d) The
corresponding light curve. In (c) and (d), the light-green and red dotted
curves show the contributions from the Compton scattering of blackbody photons
and the SSC, respectively.

We assume $R_{*}=10^{11}$ cm and $T_{*}=10^{5}$ K (Koljonen & Maccarone,
2017). Notably, that study also reveals that the photosphere extends
significantly beyond the WR star, implying that the $\gamma$-ray emitter might
reside within an optically thick wind. This would then influence the blackbody
Compton emission, as well as the $\gamma$–$\gamma$ absorption, adding
complexity to our modeling approach. To address this, we make a simplifying
assumption that the jet partially clears its surrounding environment, pushing
the optically thick wind to further distances. Thus, we neglect this
complication, which was also done in the previous studies (Dubus et al. 2010,
Z18). We first include only Compton scattering of blackbody photons. Our
results are shown in Figures 4(a,b) and Table 1. Confirming the previous
results, e.g., Z18, we find the steady state electron distribution to be
steep, with $p\approx 4.2\pm 0.1$, with a large low-energy cutoff,
$\gamma_{\rm min}\approx(2.9\pm 0.2)\times 10^{3}$, and the emission height of
$h\approx(2\pm 0.1)a$. We find that the high-energy cutoff is not required,
and set it to a high value of $\gamma_{\rm max}=10^{6}$. This implies that the
decline due to KN effects alone is sufficient to explain the steepening at
high energies in the observed spectrum. We show the folded light curve fitted
by our model vs. $\varphi$, i.e., the phase with respect to the actual
superior conjunction, assuming the best fit offset of A22. We see that the
peak of the modulation is still significantly before the superior conjunction,
implying that the jet is inclined, see Table 1. The value of the $\varphi_{\rm
j}>180\arcdeg$ means that the jet is slightly bent in the direction of the
rotation around the superior conjunction, which causes the peak of the
scattered flux to occur before it. The jet viewing angle is low, $i_{\rm
j}\approx 5\arcdeg^{+12}_{-3}$, which is in agreement with radio constraints
(Mioduszewski et al., 2001; Miller-Jones et al., 2004). The projection of the
orbital axis on the sky shows a large offset with respect to that of the jet,
$\Delta\lambda\approx 228\arcdeg^{+118}_{-35}$.

An earlier study by Zdziarski et al. (2012c) found the magnetic field strength
limited to $B\lesssim 100$ G in the $\gamma$-ray emitting region of the jet,
obtained relying only on the spectral measurements. Here, we impose stringent
constraints on the magnetic field strength employing the current $\gamma$-ray
data with the spectrum and the light curve, an upper limit for the X-ray flux
at 100 keV, and the millimeter measurements.

We have four independent constraints from the data that allow us to infer the
maximum allowed magnetic field strength. The first follows from the
$\gamma$-ray modulation light curve, for which the SSC flux cannot exceed its
minimum value. Then, the SSC contribution cannot overproduce the emission at
$\lesssim 0.1$ GeV. Depending on the magnetic field, this low-energy emission
can become quite strong given the presence of cooled low-energy electrons
below $\gamma_{\rm min}$ (see Section 3.2). The remaining two constraints are
described in Section 3.4, and are given by the observed 100 keV and 225 GHz
fluxes, which cannot be exceeded by the model spectrum.

Following this approach, we derive an upper limit for $B$ using the jet
opening angle $\alpha_{\rm j}=5\arcdeg$, a value found on the scale of
$\sim\\!10^{15}$ cm ($5.0\pm 0.5\arcdeg$; Miller-Jones et al. 2006).
Specifically, we obtain $B_{\rm max}\approx 34$ G. This magnetic field
strength yields an SSC contribution saturating the minimum of the light curve,
while we find that the X-ray and synchrotron emissions remain significantly
below the observational constraints. Figures 4(c,d) illustrate the fits to the
$\gamma$-ray data. The resulting fit parameters are given in Table 1, and are
similar to those obtained assuming zero magnetic field.

The minimum length of the emitting region can be set by $h_{\rm cool}\sim
t_{\rm cool}(\gamma_{\rm min})\beta_{\rm j}c$, where the cooling time, $t_{\rm
cool}(\gamma)\equiv\gamma/|\dot{\gamma}_{\rm cool}|$, depends on the total
cooling rate, $\dot{\gamma}_{\rm cool}$. We calculate $\dot{\gamma}_{\rm
cool}$ due to the stellar Compton process, $\dot{\gamma}_{\rm cool,SC}$, using
eq. 16 in Zdziarski et al. (2014), while the cooling rate due to the SSC (for
which the KN effects are negligible) and the synchrotron processes are
evaluated as $|\dot{\gamma}_{\rm cool,SSC/syn}|=(4\sigma_{\rm
T}\gamma^{2}U_{\mathrm{syn}/B})/(3m_{\rm e}c)$, where $U_{\mathrm{syn}/B}$ is
the energy density of the target synchrotron emission and of the magnetic
field, respectively, and $m_{\rm e}$ is the electron mass. For the model with
the maximum magnetic field we obtain $t_{\rm cool}(\gamma_{\rm min})\approx
7.7$ s. This implies $h_{\rm cool}\approx 1.6\times 10^{11}\,{\rm cm}\approx
0.25h<\Delta h=h/3$, as required (since the dissipation region is likely to be
extended). Next, we evaluate the angle-integrated jet bolometric luminosity
for the same case, obtaining $L_{\rm j}\approx 2.4\times 10^{37}$ erg s-1. The
ratio $E_{\rm e}/L_{\rm j}$ gives the average energy loss time scale of all
electrons. We find that this time scale is similar to the cooling one at
$\gamma_{\rm min}$.

###  4.2 Caveats for Model 1

Figure 5: Comparison of the orbital modulation profiles during different time
intervals. We see the profiles remain approximately constant during $\sim$5000
days of the LAT monitoring. Here, the horizontal axis gives the phase, $\phi$,
according to the ephemeris of Antokhin & Cherepashchuk (2019).  Figure 6: The
predicted angular difference between the projections on the sky of the jet and
the orbital axis, $\Delta\lambda$, see Equations (6) and (12), for the case of
a jet aligned with the spin axis of the compact object but misaligned with
respect to the orbital axis. The jet azimuthal angle, $\varphi_{\rm j}$, is
predicted to change by $360\arcdeg$ during the precession period of 50 yr,
which is clearly not seen. We assumed the best-fit values of $\theta_{\rm j}$,
$\phi_{\rm j}$ and $i$ for the case with the maximum magnetic field$B_{\rm
max}=34$ G. Currently, that model gives $\varphi_{\rm j}\approx 188\arcdeg$,
$\Delta\lambda\approx 261\arcdeg$, see Table 1.

If the jet axis is linked to the rotation axis of the compact object, prograde
geodetic precession is predicted. The De Sitter precession period for either a
BH or a neutron star is given by (Barker & O’Connell, 1975; Apostolatos et
al., 1994),

$P_{\rm prec}=\frac{c^{2}(M_{*}+M_{\rm c})^{4/3}P^{5/3}}{(2\pi
G)^{2/3}(2+3M_{*}/2M_{\rm c})M_{*}M_{\rm c}},$ (12)

where the Kepler law was used, $M_{\rm c}$ includes the contribution from the
rotational energy, and $G$ is the gravitational constant. For the best-fit
masses of A22, $P_{\rm prec}\approx 50\,{\rm yr}\approx 18000$ d, relatively
weakly depending on them. During $P_{\rm prec}$, the jet azimuth,
$\varphi_{\rm j}$, changes by $360\arcdeg$. The observations of $\gamma$-rays
by the LAT covered already $\approx 5000$ d (see Figure 1), so we would expect
substantial changes in the modulation profile. However, we see no evidence for
such changes, as shown in Figure 5.

Another argument against the presence of the precession is from the beat
between the orbital and precession periods. The two angular frequencies add,
leading to a shift of the frequency of the observed orbital modulation of the
jet $\gamma$-ray emission. This shift is $\approx$0.2 s, while the
$\gamma$-ray period measured by us equals the orbital one with the accuracy of
$<0.02$ s, see Section 2.

A consequence of the geodetic precession is also a full rotation of the jet
position angle with $P_{\rm prec}$. The predicted values of $\Delta\lambda$
are shown in Figure 6. However, we see no evidence for such changes in radio
observations from 1985 to 2016 (Molnar et al., 1988; Schalinski et al., 1995;
Mioduszewski et al., 2001; Miller-Jones et al., 2004; Egron et al., 2017). The
position angle of the approaching jet remains relatively constant at close to
$180\arcdeg$ with respect to the north.

Furthermore, this model predicts no Doppler orbital modulation of the
synchrotron emission, since it is coming from a jet with a fixed orientation.
This may disagree with observations if the preliminary results showing strong
orbital modulation at 225 GHz similar to those of the $\gamma$-rays
(McCollough, 2023) are confirmed.

The jet could possibly be misaligned due to some asymmetry in the binary,
e.g., the presence of a bow shock (A22), instead of being in the direction of
the (misaligned) BH spin. However, strong X-ray linear polarization has been
discovered from this system (Veledina et al., 2024). The polarization degree
of 21% and the polarization angle perpendicular to the position angle of the
radio jet demonstrate we do not see the primary X-ray source directly but
instead we see only the X-rays scattered from the inner back wall of an
accretion funnel (expected for super-critical accretion). This also shows that
the funnel, located close to the BH, is aligned with the jet. Thus, a
constraint on this origin of the misalignment is that it has to already start
on that scale.

##  5 Results for Model 2: orbital precession

###  5.1 Fit results

Table 2: Fit results for the phase-averaged spectrum and the light curve for a bent and precessing jet (Model 2, pure SSC) for the maximum value of the magnetic field, $B=100$ G.  $B$(G) | $\alpha_{\rm j}(\arcdeg)$ | $h/a$ | $i(\arcdeg)$ | $\theta_{\rm j}(\arcdeg)$ | $\varphi_{\rm o}(\arcdeg)$ | $\beta_{\rm j}$ | $\chi_{\nu}^{2}$ | $p$ |  $\gamma_{\rm min}$ ($\times 10^{3}$) | $\gamma_{\rm max}$ |  $E_{\rm e}\,(10^{38}$erg) | $\chi_{\nu}^{2}$ | $\gamma_{\rm br}$ | $F_{\rm X}$ | $\beta_{\rm eq}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
100 | 5(f) | 11(f) | 30(f) | $41\pm 2$ | $142\pm 4$ | $0.46\pm 0.02$ | 10.3/7 | $3.97\pm 0.06$ | $1.1\pm 0.2$ |  $10^{6}$(f) | $24\pm 3$ | 8.1/7 | 95 | 2.8 | 30.7

\begin{overpic}[height=134.42113pt]{fit_SED_CygX3_JetBend_pure_SSC_1Ledd.jpg}
\put(0.08,0.275){\large(a)} \end{overpic}

\begin{overpic}[height=134.42113pt]{fit_LC_CygX3_JetBend_pure_SSC_1Ledd.jpg}
\put(0.2,0.27){\large(b)} \end{overpic}

\begin{overpic}[height=130.08731pt]{Synchr_SED_CygX3_JetBend_pure_SSC_1Ledd.jpg}
\put(0.08,0.264){\large(c)} \end{overpic}

\begin{overpic}[height=130.08731pt]{fit_SED_CygX3_JetBend_TWO_HALVES_pure_SSC_1Ledd.jpg}
\put(0.08,0.264){\large(d)} \end{overpic}

Figure 7: Best-fit solutions for the case of the jet precessing at the orbital
period, for the derived upper limit for the magnetic field, $B_{\rm max}=100$
G. (a) The phase-averaged spectrum and (b) the 0.1–100 GeV light curve. The
blue solid curves represent the best-fit model, the dashed histograms show the
model predictions for the binned fluxes, and the jet and counterjet
contributions are depicted with the dashed cyan and dash-dotted grey curves,
respectively. (c) The predicted contribution from the synchrotron emission
along with the average flux from the SMA, shown by the black triangle. The
maximum magnetic field corresponds to the model 225 GHz flux equal to the
observed one. The jet and counterjet contributions are shown with dashed
turquoise and dash-dotted green curves, respectively. (d) Comparison of the
Fermi-LAT data and the phase-averaged models for the data split into two phase
parts corresponding to the flux higher and lower than $6\times 10^{-10}$ erg
cm-2 s-1 level (shown in blue and green, respectively).

In this model, the synchrotron and SSC emissions are orbitally modulated due
to the varying Doppler factor associated with the variable viewing angle of
the jet. The emission peaks at the lowest jet viewing angle, which occurs at
$\varphi=\pi+\varphi_{\rm o}$ (at the inferior conjunction for $\varphi_{\rm
o}=0$). Since the Doppler boosting is independent of the emission height, the
fitted SSC spectrum and modulation profile do not depend on $h$. However, the
synchrotron flux does depend on it, see Equation (7). Since both the SSC and
synchrotron are similarly modulated due to this process, this model could
explain the preliminary report of a 225 GHz orbital modulation (McCollough,
2023).

When permitting all parameters in the light curve fit to vary independently,
we observe significant uncertainties in their inferred values, indicating the
presence of degeneracies among various parameters. Specifically, the same jet
viewing angle $i_{\rm j}$ and hence the same Doppler factor modulation profile
can be achieved with different combinations of $\theta_{\rm j}$ and $i$,
implying a correlation between these parameters. To address this, we opt to
fix $i=30\arcdeg$ (A22).

As previously, we use $\alpha_{\rm j}=5\arcdeg$ (Miller-Jones et al., 2006).
We search for solutions that also reproduce the average flux measured at 225
GHz (McCollough, 2023). We find that configurations assuming a predicted
position of the jet–wind recollimation shock of $h\sim a$ (e.g., Yoon et al.
2016) significantly underpredict the 225 GHz flux, regardless of the magnetic
field strength. Consequently, we consider $h/a\gg 1$, which increases the
synchrotron flux, see Equation (7). We find $h=11a$ and $B=100$ G as the
(approximate) minimum distance of the emitting region and the maximum magnetic
field strength, respectively, yielding the observed flux at 225 GHz. The fit
parameters for this $B$ are given in Table 2 and the spectrum and the light
curve are shown in Figures 7(a,b). The corresponding model prediction for the
phase-averaged synchrotron emission is shown in Figure 7(c).

We then compare our solution to the dependence of the total magnetic field
strength in the comoving frame on the distance, $B(h)$, of the analytical jet
model of Zdziarski et al. (2022). For it, we assume the jet power of $P_{\rm
j}=9\times 10^{38}$ erg s-1, the terminal magnetization parameter of
$\sigma_{\rm min}=0.008$, where $\sigma\equiv(B^{2}/4\pi)/(\rho c^{2})$,
$\rho$ is the comoving mass density, and the spin parameter of $a_{*}=1$ and
other parameters defined in that paper of $\ell=0.5$, $a_{r}=1$,
$a_{\Gamma}=0.2$, $q_{2}=1.9$ where assumed. They yield the jet half-opening
angle of $\alpha_{\rm j}\approx 5\arcdeg$ and $B\approx 100$ G (dominated by
the toroidal component) at $h/a=11$.

In addition, we generate model predictions for the average spectra
corresponding to the Fermi-LAT high and low-flux states, delineated by the
0.1–100 GeV light curve flux being above and below the threshold of $6\times
10^{-10}$ erg cm-2 s-1, respectively. This threshold approximately corresponds
to the average between the maximum and minimum fluxes observed in the light
curve. This partitioning yields two distinct phase intervals: (1)
$0\leq\varphi/(2\pi)<0.08\ \cup 0.7\leq\varphi/(2\pi)<1$ where the 0.1 - 100
GeV flux is above the separation threshold, and (2)
$0.08\leq\varphi/(2\pi)<0.7$ where the flux is below the threshold. The
comparison of the two spectra with the model predictions is shown in Figure
7(d). We see no noticeable spectral changes between the two substates, which
are also well described by our model.

Furthermore, we estimate upper limits on $R_{*}$ and $T_{*}$ by including the
blackbody Compton component into the model. Based on observations and stellar
modeling, Koljonen & Maccarone (2017) found $T_{*}\approx 0.8$–$1.0\times
10^{5}$ K at the range of $R_{*}\approx 1.5$–$1.0\times 10^{11}$ cm,
corresponding to $R_{*}T_{*}^{2}\approx 1.0\times 10^{21}$ cm K2.
Interestingly, we find a non-monotonic behavior in the $\chi^{2}$ of the light
curve fit as a function of $R_{*}T_{*}^{2}$; it is initially decreasing, then
increasing again. The best fit is achieved for a value of
$R_{*}T_{*}^{2}=1.6\times 10^{21}$ cm K2 with $\chi_{\nu}^{2}=7.8/7$
($\Delta\chi^{2}=2.5$ with respect to the pure SSC case). However, the
$\chi^{2}$ of the spectral fit increases monotonically with increasing
$R_{*}T_{*}^{2}$, with an increase by $\Delta\chi^{2}=2.7$ as compared to the
pure SSC fit being achieved coincidentally at the same value of
$R_{*}T_{*}^{2}$ ($\chi_{\nu}^{2}=10.8/7$). Given the best fit of the light
curve and a still satisfactory SED fit, this value of $R_{*}T_{*}^{2}$ can be
therefore considered as the most optimal one. The best-fit range for
$R_{*}T_{*}^{2}$ for this most optimal solution found by imposing
$\Delta\chi^{2}\leq 2.7$ for the spectral fit (as it degrades much quicker),
is $R_{*}T_{*}^{2}=(0$–$1.7)\times 10^{21}$ cm K2. The higher $h$, the higher
the allowed range of $R_{*}T_{*}^{2}$. In addition, within this best-fit
range, the model-predicted X-ray flux at 100 keV shows a slow decrease with an
increasing $R_{*}T_{*}^{2}$, thus having consistently lower values compared to
the pure SSC case (due to adjusted $E_{\rm e}$ in the fit). We find the X-ray
flux at 100 keV for $R_{*}T_{*}^{2}=1.6\times 10^{21}$ cm K2 of $\approx
1.4\times 10^{-10}$ erg cm-2 s-1, which is below the observational constraint.

###  5.2 Jet parameters

The model of Zdziarski et al. (2022) gives the magnetic flux threading the
black hole as a function of $P_{\rm j}$ in their equation (17), which is based
on numerical results of Tchekhovskoy et al. (2009). It yields $\Phi_{B}\approx
5\times 10^{21}$ G cm2. On the other hand, since Cyg X-3 appears to be a
super-Eddington accretor (Veledina et al., 2024), $\dot{M}_{\rm accr}\gtrsim
2\times 10^{19}$ g s-1 (using the Eddington luminosity for He and the
accretion efficiency of 0.1). The magnetic flux at the limit of the
magnetically arrested disk (MAD; Bisnovatyi-Kogan & Ruzmaikin 1974; Narayan et
al. 2003) is $\Phi_{\rm MAD}\approx 50(\dot{M}_{\rm accr}c)^{1/2}r_{\rm g}$
(Tchekhovskoy et al., 2011), implying $\Phi_{\rm MAD}\gtrsim 4\times 10^{22}$
G cm2. Thus, the jet has the power at least an order of magnitude below the
MAD limit.

We then find the electron cooling time at $\gamma_{\rm min}$ of $t_{\rm
cool}(\gamma{\rm min})\approx 23$ s, translating to the cooling length of
$\approx 3.2\times 10^{11}$ cm $\approx 0.1h<\Delta h=h/3$, consistent with
our assumption of the value of $\Delta h$. We find the total radiative power
of $L_{\rm j}\approx 1.2\times 10^{38}$ erg s-1. The average energy loss time
of all electrons, $E_{\rm e}/L_{\rm j}$ (Table 2) is then in good agreement
with the cooling time scale.

Also, we calculate the Thomson optical depth, $\tau=n_{\rm e}\sigma_{\rm T}R$,
where $n_{\rm e}$ is the number density of electrons from $\gamma=1$ to
$\gamma_{\rm max}$. We find a low value of $\tau\approx 3.7\times 10^{-5}$.

We then calculate the power in both the jet and counterjet. We use eqs.
(19–21) in Zdziarski (2014). We assume a homogeneous emitting region (without
clumping), pure electron-ion plasma (fully ionized He) without positrons, and
a toroidal magnetic field. We first compute the jet power in electrons
(excluding their rest energy). We obtain $P_{\rm e}\approx 1.2\times 10^{38}$
erg s-1. Next, we evaluate the total jet power in (cold) iond, as $P_{\rm
i}\approx 5.1\times 10^{38}$ erg s-1. The magnetic power is much lower,
$P_{B}\approx 3.8\times 10^{36}$ erg s-1. The total jet power is then $\approx
6.3\times 10^{38}$ erg s-1, only slightly lower than that assumed earlier
using the model of Zdziarski et al. (2022).

Next, we compute the equipartition parameter, characterizing the degree of
equipartition between particles and magnetic field and defined in terms of the
energy densities (excluding the rest energy of particles), $\beta_{\rm
eq}=u_{\rm e}/(B^{2}/8\pi)$, where $u_{\rm e}=E_{\rm e}/(\pi R^{2}\Delta h)$
is the electron energy density. We find $\beta_{\rm eq}\approx 31$. Similarly,
the magnetization parameter (defined with respect to the rest mass density,
see above) is very low, $\sigma\sim 0.001$. This indicates the jet is strongly
dominated by the kinetic energy of particles. However, this value of $\sigma$
is lower than that assumed in the model of Zdziarski et al. (2022),
$\sigma_{\rm min}=0.008$, which is due to the magnetic power in our fit lower
than that following from that model.

In addition, we also compute the mass flow rate in ions in the jet (equation
17 in Zdziarski 2014), obtaining $\dot{M}_{\rm j}\approx 4.6\times 10^{18}$ g
s-1. We compare this rate to the super-Eddington accretion rate mentioned
earlier, $\dot{M}_{\rm accr}\gtrsim 2\times 10^{19}$ g s-1, and the derived
jet powers to $\dot{M}_{\rm accr}c^{2}\approx 1.8\times 10^{40}$ erg s-1. We
find that the jet is powered by a rather sizeable fraction ($\approx 23$%) of
accretted matter, while only a very small fraction ($\approx 3$%) of the
energy released from accretion ($\dot{M}_{\rm acc}c^{2}$) is supplied to the
particle energy, and an almost negligible fraction is converted into the jet
magnetic energy.

Next, we calculate the jet bending angle, $\Phi$, owing to the wind impact. In
the limit of $\Phi\ll 1$ rad, we have

$\Phi\approx\frac{\alpha_{\rm j}\dot{M}_{\rm w}v_{\rm w}(\Gamma_{\rm
j}-1)c}{4\pi\beta_{\rm j}\Gamma_{\rm j}P_{\rm j}},$ (13)

as given by equation (7) of Bosch-Ramon & Barkov (2016) (see their equation 8
for an expression valid at any $\Phi$). Here $\dot{M}_{\rm w}$ and $v_{\rm v}$
are the wind mass loss rate and velocity, respectively, and $\Gamma_{\rm j}$
is the jet Lorentz factor. We assume $\dot{M}_{\rm w}\approx
10^{-5}{M}_{\sun}$ yr-1 (A22) and $v_{\rm w}\approx 1.5\times 10^{8}$ cm s-1
(van Kerkwijk et al., 1996). However, since the jet likely undergoes a
decelaration over the $\gamma$-ray emission region and possibly upstream of
it, we use $\Gamma_{\rm j}=2$. This yields a low angle of $\Phi\approx
1\arcdeg$.

###  5.3 Caveats for Model 2

The main problem for this model is the value of the jet wind-induced bending
angle, $\Phi\approx 1\arcdeg$, being much lower than the fitted value of
$\theta_{\rm j}\approx 41\arcdeg$. If the (unknown) value of $\Gamma_{\rm j}$
is larger, we can increase $\Phi$, but only to $\approx\\!1.8\arcdeg$ for
$\Gamma_{\rm j}\gg 1$. Then, the jet half-opening angle in the bending region
can be larger, as seen in the hard state (Yang et al., 2023). Then, the jet
power could be lower depending on the poorly known jet parameters. In
particular, if the jet composition consists mostly of e± pairs, $P_{\rm i}$
can be negligible compared to $P_{\rm e}$, which would then increase $\Phi$ by
a factor of $\approx$4\. Furthermore, values of $\dot{M}_{\rm w}$
significantly higher than $10^{-5}{M}_{\sun}$ yr-1 were obtained by van
Kerkwijk (1993) and Ogley et al. (2001). In principle, $\Phi\approx 40\arcdeg$
could be reached if $\alpha_{\rm j}$, $\dot{M}_{\rm w}$ and $\Gamma_{\rm j}$
are higher than the adopted values, and $P_{\rm j}$ is lower.

In addition, if $\Phi<\alpha_{\rm j}=5\arcdeg$, there would be only a minor
modification of the jet propagation due to the thrust of the stellar wind,
which would only induce a weak recollimation shock or a sound wave in the jet
(Bosch-Ramon & Barkov, 2016). However, we argued above that $\Phi$ could be
larger, and it appears it could be larger than $\alpha_{\rm j}$, if, in
particular, the jet consists of a significant number of e± pairs.

A major consideration for this model is the required large phase offset,
$\varphi_{\rm o}\approx 142\arcdeg\pm 4\arcdeg$. For our assumed masses, the
BH velocity around the center of mass is $v_{\rm BH}\approx 6\times 10^{7}$ cm
s-1. The average vertical velocity of the jet between the origin and the
height of the emission required (via the binary rotation) to reach
$\varphi_{\rm o}$ is then $h/(P\varphi_{\rm o}/2\pi)$. Obviously, it is rather
low, e.g., $\approx 4.3\times 10^{8}$ cm s${}^{-1}\sim 0.01c$ for $h=11a$.
This may still possible be if the jet is composed of a slow sheath, formed by
a disk outflow, and a fast spine powered by the spin of the rotating BH.

Another issue is the large fitted polar angle of the orbital precession,
$\theta_{\rm j}=41\arcdeg\pm 2\arcdeg$ in the $\gamma$-ray emission region.
The jet needs to recollimate above that region in order to satisfy the
observational constraint on the semi-opening angle, $\approx\\!5\arcdeg$
(Miller-Jones et al., 2004, 2006). This could be achieved by reconfinement,
for instance one induced by the stellar wind (e.g., Yoon & Heinz 2015; Yoon et
al. 2016), and the symmetry and pressure gradients of the environment at the
scales of interest. In fact, the large scale jet half-opening angle can become
moderate once its interaction with the wind is not significant any more (see,
e.g., equation 19 in Bosch-Ramon & Barkov 2016). However, due to many scales
involved, estimating the jet large-scale properties would require numerical
simulations, which are out of scope of this work.

##  6 Discussion

After exhaustive research, we are still unable to choose between the two
studied models of the observed $\gamma$-ray modulation on the orbital period.
The blackbody Compton model gives better fits, but it predicts general-
relativistic jet precession on the period of 50 years, which is ruled out
observationally, see Section 4.2. The model with with a jet precessing on the
orbital period due to bending by the thrust of the stellar wind fits the data
satisfactorily, but the wind is too weak and/or the jet too strong to account
for the fitted large bending angle, see Section 5.3. This phenomenon has been
theoretically predicted (e.g., Yoon & Heinz 2015; Yoon et al. 2016; Bosch-
Ramon & Barkov 2016; Barkov & Bosch-Ramon 2022), but not yet confirmed
observationally in any other source.

An important test possibly allowing us to distinguish between the two models
would be either the confirmation or disproval of the presence of orbital
modulation with the folded light curve similar to that seen in $\gamma$-rays
at the mm-wavelength modulation. Such modulation was claimed in a preliminary
report of (McCollough, 2023). If present, it clearly cannot be explained by
the Compton anisotropy effect present in our Model 1, but it can be due to the
Doppler boosting present in our Model 2. However, we note that modulation at
mm wavelengths could also be due to free-free absorption in the stellar wind,
and the orbital phase of the maximum absorption can be shifted due to the jet
misalignment with respect to the orbital axis.

We consider then the origin of the structures observed by Mioduszewski et al.
(2001) and Miller-Jones et al. (2004), which, if caused by precession, would
have much longer periods than the binary one, $\gtrsim$60 d and $5.5\pm 0.5$
d, respectively. We speculate that these long periods are the result of
perturbations in the jet on large scales triggered by instability growth as
the jet is being perturbed at its base by the effect of the wind and the
orbital motion. This would give rise to a distinct pattern not synchronized
with the orbital cycle but rather representing a subharmonic of it.
Specifically, if the jet bending is excited with a wavelength $N$ times the
orbital period multiplied by the jet velocity, the resulting spatial scale
would be $\sim 10^{3}Na$.

A consequence of Model 2 is that the jet (on average) and the binary should
have the same inclination. The binary inclination determined by A22 is
$i=29\fdg 5\pm 1\fdg 2$. The jet inclination determined by Miller-Jones et al.
(2004) based on the proper motion of the approaching and receding components
assuming they are intrinsically symmetric (Mirabel & Rodríguez, 1994) and
$D=9$ kpc is $i_{\rm j}\approx 26\arcdeg\pm 6\arcdeg$, i.e., consistent with
the above value of $i$. On the other hand, the jet inclination based on
fitting the precessing jet model of Hjellming & Johnston (1981) to the
curvature observed in the VLBA images of Cyg X-3 is $i\approx 10\fdg 4\pm
3\fdg 8$ (Miller-Jones et al., 2004) and $i\lesssim 14\arcdeg$ (Mioduszewski
et al., 2001). Those two estimates are incompatible with the inclination of
A22, which would then favor Model 1.

In Model 2, we considered the $\gamma$-ray emission region at a relatively
large distance from the origin, $h\approx 10a$, which allowed us to obtain the
synchrotron emission strong enough to account for the observed average 225 GHz
flux. This distance is much larger than that expected for a recollimation
shock from the jet-wind interactions, which is at $h\sim(1$–$2)a$ (Yoon et
al., 2016; Bosch-Ramon & Barkov, 2016; Barkov & Bosch-Ramon, 2022). Then, the
dissipation at $h\sim 10a$ would need to be explained by another process.
Alternatively, we can relax the requirement of reproducing the 225 GHz flux by
assuming that emission originates in a different part of the jet. Then, the
SSC emission alone in that model is independent of $h$, and we could have
$h\sim(1$–$2)a$.

##  7 Conclusions

Our main conclusions are as follows.

We have used a $\gamma$-ray bright data set from the LAT observations of Cyg
X-3 greatly enlarged with respect to the previous study of Z18. We have found
that the $\gamma$-rays are modulated at the period that is compatible with the
orbital period based on the X-ray data at the accuracy of $<0.02$ s
(accounting for the secular period increase).

We have studied possible models of the orbital modulation. We have
significantly improved the model based on anisotropy of Compton scattering, by
including KN effects and the presence of a magnetic field. The model fits well
both the modulation profile and the average spectrum, and requires that the
jet is misaligned with respect to the orbital axis by $\sim\\!35\arcdeg$. If
the jet direction follows the BH spin axis, both will undergo the geodetic
precession with a period of $\sim$50 yrs. Its presence is ruled out by both
the $\gamma$-ray and radio data. If this is the correct model, the jet
misalignment is not related to the BH spin axis.

Therefore, we have proposed an alternative model using the jet bending due to
the thrust of the stellar wind from the donor. The jet is aligned with the
orbital axis on average, but it is bent to outside of the orbit, causing its
precession at the orbital period. The Doppler boosting of both the synchrotron
and SSC emission is then orbitally modulated. The modulated SSC emission fits
well the $\gamma$-ray data. The modulation of the SSC in this model does not
depend on the distance of the emission region from the center. However, the
synchrotron flux (at the SSC flux determined by the $\gamma$-ray data) depends
on it. If we postulate that the observed average mm flux is reproduced by this
model, this distance is about 10 orbital separations. The mm emission can
still be from a different jet region, and then the $\gamma$-ray emission
distance can be lower. However, a major problem for this model is that the
theoretically predicted bending angle at plausible wind and jet parameters is
much lower than that fitted to the data.

In both models, the jet velocity is low, $\beta_{\rm j}\sim 0.5$, and the jet
is weakly magnetized. The fitted field strengths are well below that predicted
by models of the magnetically arrested disk and the magnetically launched jet.

## Acknowledgements

We thank I. Antokhin, E. Egron, M. Liska and M. McCollough for discussions. We
acknowledge support from the Polish National Science Center grants
2019/35/B/ST9/03944, 2023/48/Q/ST9/00138, and from the Copernicus Academy
grant CBMK/01/24. The work of AD is supported by the South African Department
of Science and Innovation and the National Research Foundation through the
South African Gamma-Ray Astronomy Programme (SA-GAMMA).

## References

  * Antokhin & Cherepashchuk (2019) Antokhin, I. I., & Cherepashchuk, A. M. 2019, ApJ, 871, 244, doi: 10.3847/1538-4357/aafb38
  * Antokhin et al. (2022) Antokhin, I. I., Cherepashchuk, A. M., Antokhina, E. A., & Tatarnikov, A. M. 2022, ApJ, 926, 123, doi: 10.3847/1538-4357/ac4047
  * Apostolatos et al. (1994) Apostolatos, T. A., Cutler, C., Sussman, G. J., & Thorne, K. S. 1994, Phys. Rev. D, 49, 6274, doi: 10.1103/PhysRevD.49.6274
  * Bardeen & Petterson (1975) Bardeen, J. M., & Petterson, J. A. 1975, ApJ, 195, L65, doi: 10.1086/181711
  * Barker & O’Connell (1975) Barker, B. M., & O’Connell, R. F. 1975, Phys. Rev. D, 12, 329, doi: 10.1103/PhysRevD.12.329
  * Barkov & Bosch-Ramon (2022) Barkov, M. V., & Bosch-Ramon, V. 2022, MNRAS, 510, 3479, doi: 10.1093/mnras/stab3609
  * Belczyński et al. (2013) Belczyński, K., Bulik, T., Mandel, I., et al. 2013, ApJ, 764, 96, doi: 10.1088/0004-637X/764/1/96
  * Bhargava et al. (2017) Bhargava, Y., Rao, A. R., Singh, K. P., et al. 2017, ApJ, 849, 141, doi: 10.3847/1538-4357/aa8ea4
  * Bisnovatyi-Kogan & Ruzmaikin (1974) Bisnovatyi-Kogan, G. S., & Ruzmaikin, A. A. 1974, Ap&SS, 28, 45, doi: 10.1007/BF00642237
  * Blandford & Payne (1982) Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883, doi: 10.1093/mnras/199.4.883
  * Blandford & Znajek (1977) Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433, doi: 10.1093/mnras/179.3.433
  * Boettcher et al. (2012) Boettcher, M., Harris, D. E., & Krawczynski, H. 2012, Relativistic Jets from Active Galactic Nuclei (Wiley), doi: 10.1002/9783527641741
  * Bosch-Ramon & Barkov (2016) Bosch-Ramon, V., & Barkov, M. V. 2016, A&A, 590, A119, doi: 10.1051/0004-6361/201628564
  * Chiaberge & Ghisellini (1999) Chiaberge, M., & Ghisellini, G. 1999, MNRAS, 306, 551, doi: 10.1046/j.1365-8711.1999.02538.x
  * Dubus et al. (2010) Dubus, G., Cerutti, B., & Henri, G. 2010, MNRAS, 404, L55, doi: 10.1111/j.1745-3933.2010.00834.x
  * Egron et al. (2017) Egron, E., Pellizzoni, A., Giroletti, M., et al. 2017, MNRAS, 471, 2703, doi: 10.1093/mnras/stx1730
  * Fermi LAT Collaboration et al. (2009) Fermi LAT Collaboration, Abdo, A. A., Ackermann, M., et al. 2009, Science, 326, 1512, doi: 10.1126/science.1182174
  * Giacconi et al. (1967) Giacconi, R., Gorenstein, P., Gursky, H., & Waters, J. R. 1967, ApJ, 148, L119, doi: 10.1086/180028
  * Hjalmarsdotter et al. (2008) Hjalmarsdotter, L., Zdziarski, A. A., Larsson, S., et al. 2008, MNRAS, 384, 278, doi: 10.1111/j.1365-2966.2007.12688.x
  * Hjalmarsdotter et al. (2009) Hjalmarsdotter, L., Zdziarski, A. A., Szostek, A., & Hannikainen, D. C. 2009, MNRAS, 392, 251, doi: 10.1111/j.1365-2966.2008.14036.x
  * Hjellming & Johnston (1981) Hjellming, R. M., & Johnston, K. J. 1981, ApJ, 246, L141, doi: 10.1086/183571
  * Katarzyński et al. (2001) Katarzyński, K., Sol, H., & Kus, A. 2001, A&A, 367, 809, doi: 10.1051/0004-6361:20000538
  * Khangulyan et al. (2018) Khangulyan, D., Bosch-Ramon, V., & Uchiyama, Y. 2018, MNRAS, 481, 1455, doi: 10.1093/mnras/sty2356
  * Koljonen et al. (2010) Koljonen, K. I. I., Hannikainen, D. C., McCollough, M. L., Pooley, G. G., & Trushkin, S. A. 2010, MNRAS, 406, 307, doi: 10.1111/j.1365-2966.2010.16722.x
  * Koljonen & Maccarone (2017) Koljonen, K. I. I., & Maccarone, T. J. 2017, MNRAS, 472, 2181, doi: 10.1093/mnras/stx2106
  * Levine et al. (1996) Levine, A. M., Bradt, H., Cui, W., et al. 1996, ApJ, 469, L33, doi: 10.1086/310260
  * Mattox et al. (1996) Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396, doi: 10.1086/177068
  * McCollough (2023) McCollough, M. 2023, Conference presentation, https://drive.google.com/file/d/1l45DBFbzkbRmTcOY6XhIDoXGJmh5I1hW/view
  * McKinney et al. (2013) McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2013, Science, 339, 49, doi: 10.1126/science.1230811
  * Miller-Jones et al. (2004) Miller-Jones, J. C. A., Blundell, K. M., Rupen, M. P., et al. 2004, ApJ, 600, 368, doi: 10.1086/379706
  * Miller-Jones et al. (2006) Miller-Jones, J. C. A., Fender, R. P., & Nakar, E. 2006, MNRAS, 367, 1432, doi: 10.1111/j.1365-2966.2006.10092.x
  * Mioduszewski et al. (2001) Mioduszewski, A. J., Rupen, M. P., Hjellming, R. M., Pooley, G. G., & Waltman, E. B. 2001, ApJ, 553, 766, doi: 10.1086/320965
  * Mirabel & Rodríguez (1994) Mirabel, I. F., & Rodríguez, L. F. 1994, Nature, 371, 46, doi: 10.1038/371046a0
  * Molina & Bosch-Ramon (2018) Molina, E., & Bosch-Ramon, V. 2018, A&A, 618, A146, doi: 10.1051/0004-6361/201833681
  * Molina et al. (2019) Molina, E., del Palacio, S., & Bosch-Ramon, V. 2019, A&A, 629, A129, doi: 10.1051/0004-6361/201935960
  * Molnar et al. (1988) Molnar, L. A., Reid, M. J., & Grindlay, J. E. 1988, ApJ, 331, 494, doi: 10.1086/166575
  * Mori et al. (1997) Mori, M., Bertsch, D. L., Dingus, B. L., et al. 1997, ApJ, 476, 842, doi: 10.1086/303667
  * Narayan et al. (2003) Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69, doi: 10.1093/pasj/55.6.L69
  * Ogley et al. (2001) Ogley, R. N., Bell Burnell, S. J., & Fender, R. P. 2001, MNRAS, 322, 177, doi: 10.1046/j.1365-8711.2001.04129.x
  * Reid & Miller-Jones (2023) Reid, M. J., & Miller-Jones, J. C. A. 2023, ApJ, 959, 85, doi: 10.3847/1538-4357/acfe0c
  * Schalinski et al. (1995) Schalinski, C. J., Johnston, K. J., Witzel, A., et al. 1995, ApJ, 447, 752, doi: 10.1086/175914
  * Szostek & Zdziarski (2008) Szostek, A., & Zdziarski, A. A. 2008, MNRAS, 386, 593, doi: 10.1111/j.1365-2966.2008.13073.x
  * Szostek et al. (2008) Szostek, A., Zdziarski, A. A., & McCollough, M. L. 2008, MNRAS, 388, 1001, doi: 10.1111/j.1365-2966.2008.13479.x
  * Tavani et al. (2009) Tavani, M., Bulgarelli, A., Piano, G., et al. 2009, Nature, 462, 620, doi: 10.1038/nature08578
  * Tchekhovskoy et al. (2009) Tchekhovskoy, A., McKinney, J. C., & Narayan, R. 2009, ApJ, 699, 1789, doi: 10.1088/0004-637X/699/2/1789
  * Tchekhovskoy et al. (2011) Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79, doi: 10.1111/j.1745-3933.2011.01147.x
  * van der Klis & Bonnet-Bidaud (1989) van der Klis, M., & Bonnet-Bidaud, J. M. 1989, A&A, 214, 203
  * van Kerkwijk (1993) van Kerkwijk, M. H. 1993, A&A, 276, L9
  * van Kerkwijk et al. (1996) van Kerkwijk, M. H., Geballe, T. R., King, D. L., van der Klis, M., & van Paradijs, J. 1996, A&A, 314, 521, doi: 10.48550/arXiv.astro-ph/9604100
  * van Kerkwijk et al. (1992) van Kerkwijk, M. H., Charles, P. A., Geballe, T. R., et al. 1992, Nature, 355, 703, doi: 10.1038/355703a0
  * Veledina et al. (2024) Veledina, A., Muleri, F., Poutanen, J., et al. 2024, Nature Astronomy, in press, arXiv:2303.01174, doi: 10.48550/arXiv.2303.01174
  * Yang et al. (2023) Yang, J., García, F., del Palacio, S., et al. 2023, MNRAS, 526, L1, doi: 10.1093/mnrasl/slad111
  * Yoon & Heinz (2015) Yoon, D., & Heinz, S. 2015, ApJ, 801, 55, doi: 10.1088/0004-637X/801/1/55
  * Yoon et al. (2016) Yoon, D., Zdziarski, A. A., & Heinz, S. 2016, MNRAS, 456, 3638, doi: 10.1093/mnras/stv2954
  * Zdziarski (2014) Zdziarski, A. A. 2014, MNRAS, 445, 1321, doi: 10.1093/mnras/stu1835
  * Zdziarski et al. (2012a) Zdziarski, A. A., Lubiński, P., & Sikora, M. 2012a, MNRAS, 423, 663, doi: 10.1111/j.1365-2966.2012.20903.x
  * Zdziarski et al. (2012b) Zdziarski, A. A., Maitra, C., Frankowski, A., Skinner, G. K., & Misra, R. 2012b, MNRAS, 426, 1031, doi: 10.1111/j.1365-2966.2012.21635.x
  * Zdziarski et al. (2012c) Zdziarski, A. A., Sikora, M., Dubus, G., et al. 2012c, MNRAS, 421, 2956, doi: 10.1111/j.1365-2966.2012.20519.x
  * Zdziarski et al. (2014) Zdziarski, A. A., Stawarz, Ł., Pjanka, P., & Sikora, M. 2014, MNRAS, 440, 2238, doi: 10.1093/mnras/stu420
  * Zdziarski et al. (2022) Zdziarski, A. A., Stawarz, Ł., Sikora, M., & Nalewajko, K. 2022, MNRAS, 515, L17, doi: 10.1093/mnrasl/slac060
  * Zdziarski et al. (2018) Zdziarski, A. A., Malyshev, D., Dubus, G., et al. 2018, MNRAS, 479, 4399, doi: 10.1093/mnras/sty1618