# A comprehensive study of high-energy gamma-ray and radio emission from Cyg X-3

Andrzej A. Zdziarski, Denys Malyshev, Guillaume Dubus, Guy G. Pooley, Tyrel Johnson, Barbara de Marco, Adam Frankowski, Maria Chernyakova, A.R. Rao

► To cite this version:

Andrzej A. Zdziarski, Denys Malyshev, Guillaume Dubus, Guy G. Pooley, Tyrel Johnson, et al. A comprehensive study of high-energy gamma-ray and radio emission from Cyg X-3. Monthly Notices of the Royal Astronomical Society, 2018, 479 (4), pp.4399-4415. 10.1093/mnras/sty1618. hal-01797105

HAL Id: hal-01797105

https://hal.science/hal-01797105v1

Submitted on 19 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A comprehensive study of high-energy gamma-ray and radio emission from Cyg X-3

Andrzej A. Zdziarski,<fnref n="1" />* Denys Malyshev,<fnref n="2" />* Guillaume Dubus,<fnref n="3" />* Guy G. Pooley,<fnref n="4" /> Tyrel Johnson,<fnref n="5" />† Adam Frankowski,<fnref n="1" /> Barbara De Marco,<fnref n="1" /> Maria Chernyakova<fnref n="6" />,<fnref n="7" /> and A. R. Rao<fnref n="8" />

Accepted 2018 June 14. Received 2018 June 12; in original form 2018 April 20

## ABSTRACT

We study high-energy $\gamma$-rays observed from Cyg X-3 by the Fermi Large Area Telescope and the 15-GHz emission observed by the Ryle Telescope and the Arcminute Microkelvin Imager. We measure the $\gamma$-ray spectrum averaged over strong flares much more accurately than before and find it well modelled by Compton scattering of stellar radiation by relativistic electrons with the power-law index of $\simeq 3.5$ and a low-energy cutoff at the Lorentz factor of $\sim 10^3$. We find a weaker spectrum in the soft spectral state but only upper limits in the hard and intermediate states. We measure strong orbital modulation during the flaring state, well modelled by anisotropic Compton scattering of blackbody photons from the donor by jet relativistic electrons. We discover a weaker orbital modulation of the 15-GHz radio emission, which is well modelled by free-free absorption by the stellar wind. We then study cross-correlations between radio, $\gamma$-ray, and X-ray emissions. We find the cross-correlation between the radio and $\gamma$-ray emission peaks at a lag less than 1 d, while we detect a distinct radio lag of $\sim 50$ d with respect to the soft X-rays in the soft spectral state.

**Key words:** acceleration of particles – gamma-rays: general – gamma-rays: stars – stars: jets – stars: individual: Cyg X-3 – X-rays: binaries.

## 1 INTRODUCTION

Cyg X-3, one of the first discovered X-ray binaries (Giacconi et al. 1967), is a unique and puzzling system. The nature of its compact object remains uncertain; Zdziarski, Mikołajewska & Belczyński (2013) considered the radial velocity measurements of Hanson, Still & Fender (2000) and Vilhu et al. (2009) as well as constraints from the donor mass-loss rate and the orbital-period change and obtained the compact-object mass range of $M_c \simeq 2.4_{-1.1}^{+2.1} M_\odot$. Koljonen & Maccarone (2017) did not confirm the measurements of Hanson et al. (2000) in their infrared (IR) spectroscopic measurements and noted that it was possible that the velocity amplitude of Hanson

et al. (2000) traced the motion of the stellar wind rather than of the star. Nevertheless, Koljonen & Maccarone (2017) found the most likely mass range of $M_c \lesssim 5 M_\odot$. Thus, the current constraints allow either a neutron star or a low-mass black hole (BH). The presence of a BH appears to be favoured by considering various aspects of the X-ray and radio emission (Hjalmarsdotter et al. 2008, 2009; Szostek & Zdziarski 2008; Szostek, Zdziarski & McCollough 2008; Koljonen et al. 2010, 2018). Also, Zdziarski, Misra & Gierliński (2010) have shown that the differences between the shapes of the X-ray spectra of Cyg X-3 in its hard spectral state and those of confirmed accreting BH binaries can be accounted for by Compton scattering in the strong stellar wind from the donor, which also would account for the lack of high frequencies in its power spectra (Axelsson, Larsson & Hjalmarsdotter 2009). On the other hand, Burke, Gilfanov & Sunyaev (2017) showed that weakly magnetized neutron-star X-ray binaries in the hard state have similar spectral properties to their BH counterparts but lower electron

<footnote n="1">¹Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, Bartycka 18, PL-00-716 Warszawa, Poland</footnote>

<footnote n="2">²Institut für Astronomie und Astrophysik Tübingen, Universität Tübingen, Sand 1, D-72076 Tübingen, Germany</footnote>

<footnote n="3">³UJF-Grenoble 1/CNRS-INSU, Institut de Planétologie et d'Astrophysique de Grenoble, UMR 5274, 38041 Grenoble, France</footnote>

<footnote n="4">⁴Cavendish Laboratory, J. J. Thomson Avenue, Cambridge CB3 0HE, UK</footnote>

<footnote n="5">⁵College of Science, George Mason University, Fairfax, VA 22030, USA</footnote>

<footnote n="6">⁶School of Physical Sciences, Dublin City University, Glasnevin, Dublin 9, Ireland</footnote>

<footnote n="7">⁷DIAS, Fitzwilliam Place 31, Dublin 2, Ireland</footnote>

<footnote n="8">⁸Tata Institute of Fundamental Research, Mumbai 400005, India</footnote>

<footnote>* E-mail: andrzej.zdziarski@gmail.com (AAZ);</footnote>

<footnote>denys.malyshev@astro-uni-tuebingen.de (DM);</footnote>

<footnote>Guillaume.Dubus@univ-grenoble-alpes.fr (GD)</footnote>

<footnote>† Present address: Naval Research Laboratory, Washington, DC 20375, USA</footnote>
temperatures and softer spectra (which confirms some previous studies). This may be compatible with Cyg X-3 hosting a neutron star.

Cyg X-3 is the only known binary in the Galaxy containing both a compact object and a Wolf-Rayet star (van Kerkwijk et al. 1992, 1996; van Kerkwijk 1993; Fender, Hanson & Pooley 1999). Given its very short period of $P \simeq 0.2$ d (unusual for a high-mass binary), it is a likely progenitor of a close double degenerate system, after the donor explodes as a supernova (Belczyński et al. 2013). The merger will then be associated with emission of gravitational waves, which has important implications for the detectability of similar extragalactic systems by LIGO and VIRGO (Belczyński et al. 2013).

The most recent distance estimate is the geometric one from dust-scattering haloes of $D \simeq 7.4 \pm 1.1$ kpc (McCollough, Corrales & Dunham 2016), which agrees well with the estimate of $7.2^{+0.3}_{-0.5}$ kpc of Ling, Zhang & Tang (2009) obtained with the same method. A similar distance range is also preferred based on considering constraints on the donor mass (Koljonen & Maccarone 2017). At this $D$, its absorption-corrected bolometric X-ray luminosity reaches several times $10^{38}$ erg s$^{-1}$ in its brightest (soft) state, i.e. it reaches the Eddington limit for a $4M_{\odot}$ BH and exceeds it for a neutron star (Zdziarski, Segreto & Pooley 2016b, hereafter ZSP16). Among X-ray binaries, Cyg X-3 is the brightest and most highly variable radio source (McCollough et al. 1999), also showing resolved jets (Mioduszewski et al. 2001; Miller-Jones et al. 2004; Tudose et al. 2007; Egron et al. 2017) on the size scale from a few up to several tens of mas (at 7 kpc, 25 mas corresponds to the projected distance of 1 light day). However, larger radio structures are also observed (Martí, Paredes & Peracaula 2001) on the scale of 1 arcsec (corresponding to 40 light days).

Its high-energy (HE) $\gamma$-ray emission has been discovered by the Fermi Large Area Telescope (LAT; Atwood et al. 2009) and AGILE (Tavani et al. 2009a) in the soft spectral state (Fermi-LAT Collaboration 2009, hereafter FLC09; Tavani et al. 2009b). The GeV power-law emission and its orbital modulation appear to be due to Compton up-scattering of the stellar emission from the companion WR star by relativistic electrons in the jet (Dubus, Cerutti & Henri 2010). Cyg X-3 is one of only two X-ray binaries that are certainly powered by accretion for which HE $\gamma$-ray emission has been detected at a high statistical significance; the other one being Cyg X-1 (also a high-mass X-ray binary, hereafter HMXB), where the $\gamma$-ray emission is, however, much weaker (Zanin et al. 2016; Zdziarski et al. 2017). Among low-mass X-ray binaries (LMXBs), a $\gamma$-ray flare from V404 Cyg was detected at an $\sim 4\sigma$ significance during the 2015 outburst (Loh et al. 2016a). The relatively strong radio and HE $\gamma$-ray emissions in Cyg X-3 may be due to interaction of the jet with the stellar wind, which is very dense near the compact object in this close Wolf-Rayet system, and subsequent formation of recollimation shocks (e.g. Yoon, Zdziarski & Heinz 2016 and references therein).

The presence of a powerful jet in soft states in Cyg X-3 is significantly different from the behaviour of accreting BH LMXBs, where a short-duration transient jet can appear during hard-to-soft transitions (Fender, Belloni & Gallo 2004), and the radio emission is strongly quenched in the soft state (Corbel et al. 2000). As shown by Koljonen et al. (2010), strong radio flares in Cyg X-3 occur during the transition from the softest (hypersoft) states to harder ones in the opposite direction and at much lower hardness ratios than those in LMXBs. On the other hand, accreting BH binaries in the hard state feature a steady compact jet, emitting partially self-absorbed synchrotron radio-mm-IR emission. This is then sim-

ilar to Cyg X-3, which has a hard state with relatively strong radio emission correlated with soft X-rays, whose correlation is similar to that in BH binaries (Corbel et al. 2013; ZSP16). Therefore, we can expect some HE $\gamma$-ray emission of Cyg X-3 in its hard state, similar to the case of Cyg X-1.

In this work, we obtain HE $\gamma$-ray spectra in the flaring, hard, intermediate, and soft states. Furthermore, we present and study 15-GHz monitoring data from the Ryle Telescope and the Arcminute Microkelvin Imager (AMI). The latter data cover the entire duration of the LAT observations analysed here. We measure and model the orbital modulation of both radio and HE $\gamma$-ray emissions. We then study correlations between the $\gamma$-rays and radio and X-ray emissions.

## 2 DATA

### 2.1 FERMI LAT data reduction

We have analysed the available Fermi/LAT data (MJD 54682–57982) coming from the direction of Cyg X-3 using the latest version of the Fermi Science Tools (v10r0p5) with the P8R2_CLEAN_V6 instrument response functions. We have used the standard value of the zenith angle cut of $z_{max} = 90^{\circ}$.

Similarly to FLC09, we have considered the presence of the nearby $\gamma$-ray pulsar PSR J2032+4127, located about 30' away from Cyg X-3. That pulsar appears to be a member of a highly eccentric ($e \sim 0.95$), long-period ($\sim 25–50$ yr), binary with a massive Be star (Lyne et al. 2015; Takata et al. 2017), and the orbital motion causes strong variations of the pulse period. The spin ephemeris is given by Lyne et al. (2015), and according to it, the pulse curve of PSR J2032+4127 is dominated by two strong peaks at phases 0.13–0.19 and 0.62–0.70. However, PSR J2032+4127 is now approaching the periastron, which resulted in strong departures of the ephemeris from that of Lyne et al. (2015) after 2015. Based on the LAT data, we have been unable to update the ephemeris. Therefore, we have been unable to use the approach of FLC09 of using only the off-pulse intervals. Instead, we have relied on taking into account the emission of the pulsar (along with other sources in the region) in our fitting of the LAT data (see below). We also note that the pulsar flux contribution at the position of Cyg X-3 is much lower than that of Cyg X-3 when it is in the flaring state, and also no GeV flux enhancement from PSR J2032+4127 was observed while approaching the periastron (Takata et al. 2017). Still, we have compared the results for the Cyg X-3 spectra in different states obtained using all the data and those subtracting the pulsar peaks according to Lyne et al. (2015) and found virtually no differences. Therefore, we have used the entire data in this work.

In order to take account of the broad Fermi/LAT point-spread function at energies studied by us (80 MeV–300 GeV), we consider a large, $25^{\circ} \times 25^{\circ}$, region of interest (ROI) around the Cyg X-3 position. We include in the modelling all sources within the ROI from the 4-yr Fermi catalogue (Fermi-LAT Collaboration 2015; 3FGL). We have used the standard templates for the Galactic (gll_iem_v06.fits) and extragalactic (iso_P8R2_CLEAN_V6_v06.txt) diffuse backgrounds. The catalogue sources were assumed to be described by the 3FGL spectral models with all parameters except the normalizations frozen to their catalogue values. In order to avoid possible systematic effects, we have also included into the model the 3FGL catalogue sources located up to $10^{\circ}$ beyond the ROI with all parameters fixed to their catalogue values. For Cyg X-3, we adopt the power-law model with the slope 2.7 reported previously in FLC09. The spectral analysis
has been performed with the python tools.<fnref n="1" /> The upper limits are calculated with the *IntegralUpperLimits* python module for detection significances of TS (test statistic; see Mattox et al. 1996) <4, which correspond to a 95 per cent (≈2σ) probability for the energy flux to be lower than the limit.

We then built the TS map in a 5° × 5° region around the position of Cyg X-3 in the 1–300 GeV energy band, see Fig. 1. We see a number of residuals along the Galactic plane, which we mark as n1–n6. Since almost all of them are at very low Galactic latitudes, where we expect the highest uncertainties, most of them seem to be diffuse residuals unaccounted for in the diffuse background model. The only residual that can be identified in any catalogue is n1, which was present in the first Fermi catalogue (Fermi-LAT Collaboration 2010), but it disappeared in 3FGL, and which appears to be associated with an H II region (Munar-Adrover, Paredes & Romero 2011). The map also reveals a weak point-like source at the catalogue position of Cyg X-3 with TS ≈ 45, corresponding to a ≳6σ detection significance. Hereafter, we use the High Energy Astrophysics Science Archive Research Center (HEASARC) catalogue position of Cyg X-3, (RA, Dec) = (308.107420; 40.957750). We also note that the Cyg X-3 position is actually consistent with that of J2032+4050 as given in the preliminary 8-yr Fermi/LAT source list<fnref n="2" /> (which appeared when the paper was in final stages of preparation).

The spectral analysis in 0.08–300 GeV was performed in a set of narrow energy bins. For each energy bin, we have iteratively removed weak (TS < 1) sources other than Cyg X-3 and then have redone the fit until no weak sources remain.

In addition, we also consider a lower energy band available to the LAT of 40–80 MeV (30–104 MeV accounting for the energy dispersion) for our brightest (flaring) spectrum, where we find an upper limit. For that energy range, we employ a method similar to that used in Zdziarski et al. (2017). This range is not covered by the standard templates for Galactic and isotropic diffuse emission. Therefore, we base our analysis for the Galactic background on three different templates, $S_LZ_6R_{20}T_\infty C_5$, $S_5Z_4R_{20}T_{150}C_5$, and $S_YZ_6R_{30}T_{150}C_2$, produced with the GALPROP code (Vladimirov et al. 2011). Those templates are known to describe Fermi/LAT data at higher energies reasonably well (Ackermann et al. 2012). The spectrum of the standard isotropic background model is available down to 34 MeV, which almost covers the analysed energies, and we employ a power-law continuation of that spectrum down to 30 MeV. We also use a low-value zenith angle cut of $z_{max} = 70°$, instead of the standard value of 90° (which we adopt at higher energies).

We have then performed timing analysis for Cyg X-3 in 1-d bins in a broad energy range of 0.1–10 GeV in a way similar to the above-described binned spectral analysis with iterative elimination of weak sources. Fig. 2(a) shows the 1-d bin light curve of the LAT detections. We find 486 d with the signal-to-noise ratios (SNR) >1 and 174 d with SNR >2. The detections with TS ≥16 and <16 are plotted in blue and cyan, respectively. We confirm most of the previous detections by the LAT and AGILE (FLC09; Tavani et al. 2009b; Williams et al. 2011; Bulgarelli et al. 2012; Corbel et al. 2012; Piano et al. 2012; Bodaghee et al. 2013; Cheung & Loh 2016; Loh et al. 2016b; Piano et al. 2016, 2017a,b; Loh & Corbel 2017a,b), which we show by the grey vertical lines. We also find a number of new detections. We then split the days with detections into two energy flux regions, the high-flux region which we call

the flaring state and the other, with detections at lower fluxes. The boundary between the two regions, $F_b$, is selected with the iterative procedure defined as following. At each iteration, we define the mean level $F_{low}$ and the standard deviation $\sigma_{low}$ of all detections defined as non-flaring at the previous iteration. If any, we mark all detections with a flux higher than $F_b = F_{low} + 3\sigma_{low}$ as flares and continue to the next iteration. The iterations stop when no new detections are attributed to the flaring state. Using this procedure, $F_b$ is found to equal $4.94 \times 10^{-10}$ erg cm$^{-2}$ s$^{-1}$, which is shown by the red horizontal line in Fig. 2(a). We find 49 d with the energy flux ≥ this limit, and we list them in Table 1.

We then divide the available LAT observations into the hard, intermediate, and soft states based on the daily-averaged data from All-Sky Monitor (ASM; Bradt, Rothschild & Swank 1993; Levine et al. 1996) on board *Rossi X-ray Timing Explorer*, the Monitor of All-sky X-ray Image (MAXI; Matsuoka et al. 2009) on board *International Space Station*, and the Burst Alert Telescope (BAT; Barthelmy et al. 2005; Markwardt et al. 2005; Krimm et al. 2013) on board *Swift*. The long-term light curves from those detectors are given in ZSP16, and they are updated in Fig. 3. In order to determine the states, we use a method similar to that in Zdziarski et al. (2012a) and in ZSP16, except that we use here the public 15–50 keV BAT data<fnref n="3" /> instead of custom data used in those papers. We convert the ASM and BAT count rates and the MAXI photon fluxes into energy fluxes using the scaling to the Crab, assuming its spectrum as given in ZSP16. Figs 2(b–c) shows the BAT flux versus the 3–5 keV ASM and 2–4 keV MAXI fluxes. The flux regions delineating the states are defined by the dashed lines. The 3–5 keV boundary of the hard state corresponds to the maximum flux with a positive correlation with the radio emission (ZSP16) and the hard/soft X-rays anti-correlation, and the 15–50 keV boundary approximately corresponds to the lowest fluxes of the hard state. For days without both soft and hard X-ray data, we interpolate the available X-ray data to infer the spectral state, as well as use the 15-GHz data and the radio/X-ray correlations as given in ZSP16. This yields 3074, 736, and 446 d with LAT coverage in the hard, intermediate, and soft state, respectively.

Then, we find 43, 2, and 4 d with γ-ray flares (defined as above) in the soft, intermediate, and hard state, respectively. The points corresponding to the flares for the days with both soft and hard X-ray coverage are shown in red in Figs 2(b–c). They show the four flaring days in the hard state, where one point appears on both panels, i.e. it has both simultaneous ASM/BAT and MAXI/BAT coverages. The occurrence of the intermediate state for two flares was determined by interpolating the X-ray data (see above); thus, those days do not appear in Figs 2(b–c). We see that while most of the days with strong γ-ray detections correspond to the soft state, there are still several detections during the intermediate and hard states, i.e. with low soft X-ray energy fluxes and high hard X-ray ones.

## 2.2 The radio data

We study here radio-monitoring data at 15 GHz from the Ryle Telescope, which cover MJD 49231–53905 (74181 measurements), and the AMI, MJD 54612–58055 (5125 measurements). The combined data set contains 79306 measurements. The AMI Large Array is the re-built and reconfigured Ryle Telescope. Pooley & Fender (1997) describe the normal operating mode for the Ryle telescope in the

<footnote n="1">¹fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/python_tutorial.html</footnote>

<footnote n="2">²https://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl8y/</footnote>

<footnote n="3">³http://swift.gsfc.nasa.gov/results/transients/CygX-3/</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The TS map (Galactic coordinates) at energies ≥1 GeV for the data within a 5° × 5° square around the position of Cyg X-3 with the 3FGL sources subtracted, with their positions marked by the green crosses. The apparently detected sources not present in 3FGL are marked with cyan circles and denoted n1–n6 (see Section 2.1). Cyg X-3 is marked by the dotted cyan circle in the centre.</figcaption></figure>

monitoring observations; the observing scheme for the AMI Large Array is very similar. The new correlator has a useful bandwidth of about 4 GHz (compared with 0.35 GHz for the Ryle), but the effective centre frequency is similar.

In order to establish the calibration parameters of the array, we use observations of a bright, nearby unresolved source interleaved with those of the main target source. Our primary calibrators were 3C 48 and 3C 286. This procedure resulted in variations in the flux calibration limited to <10 per cent from one day to another.

Throughout this work, we consider the radio emission as coming from the source associated with the accretion/outflow associated with the compact object, presumably the jet. Still, some contribution to that emission comes from the stellar wind. However, it has to be minor, given the very low radio flux levels the source achieves. Another possible contribution is from the stellar wind interacting with disc winds (Koljonen et al. 2018). In this work, we will not distinguish this contribution from the jet, given that we have only 15-GHz radio fluxes at our disposal.

## 3 LAT SPECTRA IN DIFFERENT STATES

We have calculated the average LAT spectra in the hard, intermediate, soft, and flaring states, as defined in Section 2.1, except that we have excluded the flaring days from the soft-state spectrum. Our results are shown in Fig. 4.

We have detected the source with the significance of ~40σ during the flaring state in the 0.08–15 GeV range. A power-law fit of the 0.08–10 GeV range gives a photon spectral index of Γ =

2.55 ± 0.05 at the normalization at 100 MeV of 1.9 ± 0.2 × 10⁻⁸ MeV⁻¹ cm⁻² s⁻¹, which corresponds to the energy flux above 100 MeV of 5.5 × 10⁻¹⁰ erg cm⁻² s⁻¹. This power-law spectrum is significantly more accurately determined and slightly harder than the Γ = 2.70 ± 0.25 of FLC09 (measured during MJD 54750–54820 and 54990–55045), and the integrated energy flux is slightly larger than their value of 4.0 ± 1.6 × 10⁻¹⁰ erg cm⁻² s⁻¹ at their best fits. Nevertheless, both Γ and the flux are consistent with the FLC09 values within the uncertainties.

However, there is a visible curvature in the present flaring-state spectrum, and we have also fitted it by a lognormal distribution in the form of $dN/dE = N_0(E/E_b)^{-2}10^{-\beta \log_{10}(E/E_b)}$, where $E_b$ is the peak of it in $EF_E$. (Note that this form is equivalent to a parabola in logarithmic coordinates, see Zdziarski et al. 2016a.) We obtain $E_b = 384 \pm 15$ MeV, $\beta = 0.351 \pm 0.005$, and $N_0 = 1.36 \pm 0.09 \times 10^{-9}$ cm⁻² s⁻¹. We find that the lognormal/logparabola fit is strongly preferable to the power-law fit, with $\Delta\chi^2 \simeq -59$ for adding one free parameter, which corresponds to a significance of $\sigma \simeq 7–8$ of the presence of a curvature. We note that the lognormal model is sharply cut off at >2 GeV and it is thus much below the last detected spectral point.

We still detect Cyg X-3 in the soft state outside the flaring days in the 0.2–15 GeV range at a lower flux than that in the flaring state. This spectrum is parallel to the flaring state one above ~1 GeV, but we see a hardening at lower energies. We do not detect the source in the intermediate state. Contrary to our original expectations (see Section 1), we have not detected Cyg X-3 in the hard state, obtaining stringent upper limits.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) The 1-d bin γ-ray light curve with detections by the LAT. Our criterion defining the flaring state for 1-d integration, $F(0.1–100\,\text{GeV}) > 4.94 \times 10^{-10}\,\text{erg cm}^{-2}\,\text{s}^{-1}$, is shown by the red line. The grey vertical lines correspond to the previously published detections. The large blue circles and small cyan circles with flux error bars correspond to the days with $\text{TS} \ge 16$ and $<16$, respectively. (b) The relationship between the daily-averaged energy fluxes in the 3–5 keV (ASM) and 15–50 keV (BAT) ranges. The hard and soft states are defined here by $F(3–5\,\text{keV}) < 0.5\,\text{keV cm}^{-2}\,\text{s}^{-1}$ (corresponding to the count rate $\lesssim 2.7\,\text{s}^{-1}$) and $F(15–50\,\text{keV}) < 1.1\,\text{keV cm}^{-2}\,\text{s}^{-1}$ (corresponding to the count rate $\lesssim 0.028\,\text{cm}^{-2}\,\text{s}^{-1}$), respectively. These boundaries are marked by the dotted and dashed lines, respectively. The intermediate state corresponds to both the 3–5 keV and 15–50 keV fluxes above the respective boundaries. Only points with statistical significance $>3\sigma$ (required for each of the fluxes) are shown. The thick red circles with error bars correspond to MJDs with LAT flares, defined as above. (c) The same except that the 2–4 keV flux from MAXI is used. The adopted hard-state condition (the dotted line) is here $F(2–4\,\text{keV}) < 0.4\,\text{keV cm}^{-2}\,\text{s}^{-1}$ (corresponding to the photon flux $\lesssim 0.14\,\text{cm}^{-2}\,\text{s}^{-1}$).</figcaption></figure>

Fig. 5 shows the flaring-state spectrum together with the average X-ray spectra of Cyg X-3 from RXTE (Szostek et al. 2008). In their classification, the average spectra correspond to five spectral

states from the hardest to the softest. We have compared our flaring-state spectrum to the models of Zdziarski et al. (2012b), in which relativistic electrons Compton-upscatter blackbody photons emitted
<table><caption>**Table 1.** The MJD days corresponding to the flaring state defined in Section 2.1 and Fig. 2(a), with $F_γ$ (0.1–100 GeV) > 4.94 × 10⁻¹⁰ erg cm⁻² s⁻¹. The γ-ray flux and the TS are averages over a given MJD, and the letters S, I, and H correspond to the soft, intermediate, and hard state, respectively.</caption><thead><tr><th>MJD</th><th>F<sub>γ</sub> [10<sup>-10</sup> erg cm<sup>-2</sup> s<sup>-1</sup>]</th><th>TS</th><th>State</th></tr></thead><tbody><tr><td>54725</td><td>5.2 ± 1.6</td><td>18.5</td><td>H</td></tr><tr><td>54780</td><td>5.7 ± 2.3</td><td>22.7</td><td>S</td></tr><tr><td>54781</td><td>5.5 ± 1.5</td><td>23.3</td><td>S</td></tr><tr><td>54786</td><td>6.4 ± 1.6</td><td>35.3</td><td>S</td></tr><tr><td>54809</td><td>5.8 ± 1.3</td><td>36.8</td><td>S</td></tr><tr><td>54810</td><td>5.7 ± 1.4</td><td>32.3</td><td>S</td></tr><tr><td>54812</td><td>9.4 ± 1.4</td><td>83.0</td><td>S</td></tr><tr><td>54814</td><td>6.8 ± 1.3</td><td>50.4</td><td>S</td></tr><tr><td>54991</td><td>5.6 ± 1.6</td><td>18.5</td><td>S</td></tr><tr><td>54995</td><td>7.0 ± 2.2</td><td>29.0</td><td>S</td></tr><tr><td>55002</td><td>8.9 ± 1.4</td><td>67.0</td><td>S</td></tr><tr><td>55003</td><td>6.8 ± 1.5</td><td>43.8</td><td>S</td></tr><tr><td>55023</td><td>6.3 ± 1.2</td><td>39.9</td><td>S</td></tr><tr><td>55032</td><td>7.0 ± 1.8</td><td>28.2</td><td>S</td></tr><tr><td>55034</td><td>12.0 ± 1.5</td><td>119.2</td><td>S</td></tr><tr><td>55035</td><td>6.4 ± 2.1</td><td>38.7</td><td>S</td></tr><tr><td>55043</td><td>7.2 ± 1.5</td><td>40.4</td><td>S</td></tr><tr><td>55328</td><td>5.5 ± 1.3</td><td>29.5</td><td>S</td></tr><tr><td>55341</td><td>5.7 ± 1.3</td><td>28.4</td><td>S</td></tr><tr><td>55342</td><td>9.2 ± 1.5</td><td>61.5</td><td>S</td></tr><tr><td>55343</td><td>7.5 ± 1.6</td><td>38.7</td><td>S</td></tr><tr><td>55526</td><td>5.7 ± 1.4</td><td>24.8</td><td>H</td></tr><tr><td>55592</td><td>6.5 ± 1.3</td><td>40.0</td><td>S</td></tr><tr><td>55596</td><td>6.5 ± 1.7</td><td>31.1</td><td>S</td></tr><tr><td>55600</td><td>5.2 ± 2.1</td><td>17.6</td><td>S</td></tr><tr><td>55604</td><td>7.8 ± 1.8</td><td>37.2</td><td>S</td></tr><tr><td>55605</td><td>6.2 ± 1.9</td><td>23.9</td><td>S</td></tr><tr><td>55642</td><td>6.5 ± 1.9</td><td>40.5</td><td>S</td></tr><tr><td>55888</td><td>5.0 ± 1.2</td><td>29.6</td><td>I</td></tr><tr><td>55921</td><td>5.2 ± 1.5</td><td>19.0</td><td>H</td></tr><tr><td>56649</td><td>6.6 ± 2.0</td><td>19.3</td><td>H</td></tr><tr><td>56766</td><td>6.1 ± 2.2</td><td>11.5</td><td>I</td></tr><tr><td>57367</td><td>5.1 ± 1.4</td><td>19.7</td><td>S</td></tr><tr><td>57402</td><td>6.4 ± 1.4</td><td>34.0</td><td>S</td></tr><tr><td>57414</td><td>8.0 ± 1.8</td><td>37.3</td><td>S</td></tr><tr><td>57621</td><td>8.2 ± 1.9</td><td>35.9</td><td>S</td></tr><tr><td>57631</td><td>5.8 ± 1.5</td><td>24.4</td><td>S</td></tr><tr><td>57646</td><td>6.1 ± 1.2</td><td>39.0</td><td>S</td></tr><tr><td>57647</td><td>6.6 ± 1.6</td><td>41.0</td><td>S</td></tr><tr><td>57649</td><td>11.7 ± 1.5</td><td>108.6</td><td>S</td></tr><tr><td>57799</td><td>6.5 ± 1.9</td><td>42.3</td><td>S</td></tr><tr><td>57805</td><td>5.5 ± 1.3</td><td>35.0</td><td>S</td></tr><tr><td>57810</td><td>6.3 ± 1.6</td><td>35.7</td><td>S</td></tr><tr><td>57816</td><td>5.7 ± 1.4</td><td>32.0</td><td>S</td></tr><tr><td>57818</td><td>6.0 ± 1.7</td><td>42.9</td><td>S</td></tr><tr><td>57825</td><td>6.6 ± 1.4</td><td>46.3</td><td>S</td></tr><tr><td>57826</td><td>6.2 ± 1.5</td><td>31.7</td><td>S</td></tr><tr><td>57839</td><td>8.2 ± 1.3</td><td>63.1</td><td>S</td></tr><tr><td>57852</td><td>6.1 ± 2.0</td><td>15.9</td><td>S</td></tr></tbody></table>

by the donor, and which take into account the full Klein–Nishina cross section. Among those, the model with the steady-state electron power-law index of $p = 3.5$, the minimum electron Lorentz factor of $\gamma_1 = 1300$, the maximum one of $\gamma_2 \to \infty$, and scattering stellar blackbody photons at the temperature of $10^5$ K, fits well the current spectrum with $\Delta\chi^2 \simeq -23$ with respect to the above power-law fit. Only the normalization has been fitted, yielding the flux at 1 GeV of $EF_E \simeq 0.0638$ keV cm⁻² s⁻¹. The Lorentz factor of $\gamma_1$ corresponds

to the minimum above which the electrons are accelerated with an index of $p_{acc} \simeq 2.5$. Below $\gamma_1$, the electrons are from cooling by the Compton and adiabatic losses and have the distribution given by equation (21) of Zdziarski et al. (2012b). The electron spectral index is harder than that corresponding to Compton scattering in the Thomson limit, $p = 2\Gamma - 1 = 4.1$ because of the Klein–Nishina decline of the Compton cross section, which softens the spectrum. The model satisfies the constraint obtained by Zdziarski et al. (2012a) that the contribution of the jet emission at ~100 keV is minor (based on the pattern of the orbital modulation found at 50–100 keV). We confirm the result of Zdziarski et al. (2012b) that the magnetic field strength in the γ-ray-emitting region is relatively weak, $B \lesssim 100$ G.

Fig. 5 also shows the hard-state spectrum upper limits. We can see that they are quite stringent, implying any jet emission in that state to have an $EF_E$ spectrum at a level a few thousand times below the peak of the hard-state X-ray spectrum. We note that the hard-state HE γ-ray spectrum of Cyg X-1 is actually about four orders of magnitude below the peak of the hard-state spectrum (Zanin et al. 2016; Zdziarski et al. 2017). So a γ-ray spectrum at a similar relative level could still be emitted by Cyg X-3 and remain undetectable.

# 4 ORBITAL MODULATION

## 4.1 The ephemeris

The period of Cyg X-3 is increasing. We take it into account by using a quadratic form of the ephemeris,

$$T_n = T_0 + P_0 n + c_0 n^2, \quad c_0 = P_0 \dot{P}/2, \quad P_n = P_0 + 2c_0 n, \quad (1)$$

where $T_n$ is the time of an nth occurrence of a zero orbital phase, approximately<fnref n="4" /> related to the superior conjunction, and measured from the reference time, $T_0$, $P_0$ is the period at $T_0$, $\dot{P}$ is the period derivative, and $P_n$ is the period at $T_n$. The most recent ephemeris is that of Bhargava et al. (2017),

$$T_0 = 40949.384, \quad P_0 = 0.19968476(3) \text{ d}, \\ c_0 = 5.41(2) \times 10^{-11} \text{ d}, \quad (2)$$

where hereafter the numbers in parentheses give the uncertainty of the last digit.

The above ephemeris is given in the Terrestrial Time MJD and is based on X-ray light curves taking into account the barycentric correction (Y. Bhargava, private communication). Thus, we consider the same time format and apply the barycentric correction to the light curves used for the orbital modulation. In order to determine the orbital phase of a measurement at a time $T$, we solve equation (1) for $n$ treating it as a real number and then subtract its integer part.

## 4.2 Modulation of HE γ-rays

We clearly detect the period of Cyg X-3 in HE γ-rays in the flaring state. We note that the orbital period of Cyg X-3 is increasing, which could shift and smear out the peak due to the periodicity. To account for that, we convert the observation time to that of a

<footnote n="4">⁴We note that the template of van der Klis & Bonnet-Bidaud (1989) has the minimum slightly below zero phase, which is also the case for the X-ray light curves phase-folded based on a previous ephemeris in Zdziarski et al. (2012a). Furthermore, the strongest X-ray absorption may not exactly correspond to the conjunction due to a likely asymmetry of the stellar wind in this short-period binary.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The recent light curves of Cyg X-3 normalized to the respective average over the total observation length. The blue squares, red triangles, and black crosses show the rates normalized to the averages over all available data for the MAXI (2–10 keV), BAT (15–50 keV) and AMI (15 GHz), equal to 0.37 cm⁻² s⁻¹, 0.0338 cm⁻² s⁻¹, 0.104 Jy, respectively. Only points with the significance ≥2σ are shown; for clarity of display, the error bars are not plotted. The dashed green line corresponds to the averages. The heavy magenta vertical lines correspond to the flares of HE γ-ray emission, as defined in Section 2.1. We note that the last days of the HE γ-ray and 15-GHz data analysed in this work are MJD 57982 and 58055, respectively.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** The Fermi LAT spectra and upper limits for the sum of all single MJDs with strong γ-ray detections (the flaring state, as defined in Fig. 2(a), and the soft (excluding the flaring days), intermediate, and hard states, shown by the magenta inverted triangles, red circles, green squares, and blue triangles, respectively, with associated error bars and arrows indicating upper limits. The dashed line and the shaded region show the best power-law fit to the flaring-state spectrum and its uncertainties.</figcaption></figure>

constant period by calculating *n*(*T*) for an observation time, *T*, by treating *n* as a real number and solving the binomial in equation (1). We then subtract *c*₀*n*² from the time of an observation. The results of our Lomb–Scargle analysis for the light curve corrected in this way are shown in Fig. 6. We find the period of 0.199688(4) d, which, given its standard deviation, agrees very well with *P*₀ of Cyg X-3 of equation (2). We have also found an analogous result on the direct flaring-state light curve, with the peak corresponding to the range of the orbital period during the epoch of the LAT observations, though

<figure><img src="image_7.png" /><figcaption>**Figure 5.** The Fermi LAT measurements and upper limits in the flaring state (upper γ-ray error bars) and the upper limits in the hard state (lower γ-ray error bars) compared to the average RXTE spectra of Cyg X-3 in the five spectral states of Szostek et al. (2008). The X-ray hard-state spectra are the red and blue ones, i.e. with the two top fluxes at 20 keV and bottom at 4 keV, and the remaining spectra correspond to our intermediate and soft states. The solid curve shows the Compton-scattering model spectrum from Zdziarski et al. (2012b) with the electrons with a power-law distribution with the index of *p* = 3.5 above the low-energy break at γ₁ = 1300 and extending up to γ₂ → ∞. The electrons scatter stellar blackbody photons with the temperature of 10⁵ K. The fitted normalization is larger by a factor 2.1 than that of Zdziarski et al. (2012b).</figcaption></figure>

with a lower peak power, reflecting a period change during that epoch.

We then use the ephemeris of Section 4.1 to assign the phase to each photon observed within the ROI during the flaring state and split the data over six equal phase bins. We perform the binned like-
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The Lomb–Scargle periodogram for Cyg X-3 in the flaring state in the 0.1–100 GeV range, calculated accounting for the orbital period increase, taking into account the measurement uncertainties, and normalized with the $\chi^2$ of a constant model (which results in the power within the 0–1 range). The highest peak corresponds to the orbital period of Cyg X-3.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** The observed orbital modulation in the flaring state in the 0.1–100 GeV range shown by the error bars. The zero phase is defined from X-rays (Bhargava et al. 2017) and it approximately corresponds to the superior conjunction. For clarity, hereafter two full phase ranges are shown. The blue solid line shows the best fit of the Compton-anisotropy model, and the red dashed and magenta dotted lines show the contributions of the jet and counterjet, respectively.</figcaption></figure>

lihood analysis (see Section 2.1) in each of the bins.<fnref n="5" /> The resulting energy fluxes as a function of the orbital phase are shown in Fig. 7. We have not found any statistically significant dependence on the energy range, which we looked for by using the photon energy ranges of 0.1–1, 1–10, and 10–100 GeV. We have also verified the consistency of our analysis by checking that the light curve of the nearby pulsar PSR J2032+4127 remains constant in all considered phase bins.

We then model the orbital light curve obtained by Compton anisotropy. This method utilizes two features of Compton scatter-

<figure><img src="image_10.png" /><figcaption>**Figure 8.** The geometry of Compton scattering of the blackbody photons. The axes x and y are in the binary plane, and the +z direction gives the binary axis. The +x direction gives the projection of the direction towards the observer onto the binary plane. The observer is at an angle, i, with respect to the orbital axis, $\phi$ is the orbital phase, $\phi = 0$ and $\pi$ correspond to the superior and inferior conjunction, respectively, $\theta_j$ is the inclination of the jet with respect to the binary axis, $\phi_j$ is the angle of the projection of the jet onto the binary plane with respect to x-axis, d is the distance between the stars, H is the distance of the $\gamma$-ray source from the centre of the compact object, the vectors $e_{\text{obs}}$, $e_c$, and $e_*$ point from the donor towards the observer, the centre of the compact object, and the $\gamma$-ray source, respectively, and $e_j$ points from the centre of the compact object towards the $\gamma$-ray source.</figcaption></figure>

ing of stellar emission by relativistic electrons. First, the scattering probability is maximized for head-on collisions, i.e. for electrons moving towards the star. Second, a relativistic electron emits the scattered photon predominantly along its direction of motion. Thus, most of the Compton-scattered emission is towards the star and almost no emission is directed to an observer located along the line connecting the star centre and the $\gamma$-ray source. Therefore, the observed emission is maximized when the $\gamma$-ray source is behind the star. For a jet perpendicular to the orbital plane, this would be at the superior conjunction, and the modulation would be symmetric around it. Departures from that indicate that the jet is inclined with respect to the binary axis.

We use the method of Dubus et al. (2010), including minor corrections given in Zdziarski et al. (2012b), and use the coordinate system shown in Fig. 8 (in which $\phi = 0$ corresponds to the superior conjunction). We assume the blackbody photons to be emitted by a point source with the luminosity of $L_* = 4\pi R_*^2 \sigma_B T_*^4$, where $R_*$ and $T_*$ are the stellar radius and temperature, respectively. We also assume the $\gamma$-ray source to be a point source, located at a distance, H, from the centre of the compact object. We assume the Thomson limit of Compton scattering, see equation (A9) of Zdziarski et al. (2012b). We take into account the emission of both the jet and counterjet and exclude fits to the observed modulation in which the counterjet is obscured by the star. We calculate the power injected into the non-thermal electrons in the jet+counterjet with a power-law distribution with the index of $p = 4.1$ (corresponding to the fitted power-law index in the Thomson limit, $p = 2\Gamma - 1$) and the minimum Lorentz factor of $\gamma_{\text{min}} = 10^3$. We assume fast cooling and thus that power equals the Compton-scattered lumi-

<footnote n="5">⁵We note that FLC09 performed an aperture analysis to determine the orbital modulation, which left an uncertainty about the background level, see their fig. 3B. The present method avoids this problem and relies on the standard templates, see Section 2.1.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The distributions of the mutual dependences of the found acceptable parameters (within 68 per cent confidence region), $\phi_j$ versus $\theta_j$ (top panel), and $H/d$ versus $\beta$ (bottom panel), obtained with the Compton-anisotropy model applied to the observed $\gamma$-ray orbital modulation. The colour changes from white to dark blue $\propto \sqrt{n}$, where $n$ is the number of acceptable models per unit area.</figcaption></figure>

nosity emitted by the jet in all directions. In the calculations, we assume the donor mass of $M_* = 14\text{M}_\odot$, $M_c = 4.5\text{M}_\odot$ (yielding the separation of $d = 2.65 \times 10^{11}$ cm) and the orbital inclination of $i = 31^\circ$, which correspond to the solution with the largest allowed masses in Zdziarski et al. (2013). We also assume no eccentricity, $R_* = 10^{11}$ cm, $T_* = 10^5$ K, and $D = 7$ kpc. The assumed stellar radius is less than the Roche lobe radius, $\simeq 1.16 \times 10^{11}$ cm (for the assumed masses; Eggleton 1983).

Our best-fitting solution, shown in Fig. 7, gives the jet velocity of $\beta \simeq 0.73$, the location of the $\gamma$-ray source along the jet at $H \simeq 6.1 \times 10^{11}$ cm ($\simeq 2.3d$), and the jet inclination with respect to the orbital axis of $\theta_j \simeq 37^\circ$, with an azimuthal angle $\phi_j \simeq 5^\circ$ (Fig. 8). Hence, the jet direction is $\simeq 6^\circ$ off from the line-of-sight, i.e. the jet is nearly pointed towards us. The total $\chi^2$ of the fit is 1.35 (with six orbital flux measurements and five fitted parameters, including normalization to the flux). The power injected into the non-thermal electrons (and/or $e^\pm$ pairs) is $\simeq 3.1 \times 10^{36}$ erg s$^{-1}$, and the energy content of the electrons is $\simeq 1.1 \times 10^{38}$ erg. In most of the acceptable solutions, the injected power is between $10^{36}$ and $10^{37}$ erg s$^{-1}$. Fig. 9 shows the mutual dependences between $\theta_j$ and $\phi_j$, and $\beta_j$ and $H/d$. We generally find that the jet has to be inclined with respect to the binary axis by a relatively large angle, $\theta_j \gtrsim 25^\circ$, with the $\gamma$-ray

emission zone located far from the compact object at $H \gtrsim d$. The acceptable ranges of our solutions are significantly narrower than those of Dubus et al. (2010) but still consistent with them.

We see no evidence for precession within the epoch of the studied LAT observations in either the power spectrum or by comparing the modulation shape at various epochs. It is likely that the inclined jet is aligned with the black hole spin axis up to the location of the $\gamma$-ray emission. Occasional jet precession observed in radio (Mioduszewski et al. 2001; Miller-Jones et al. 2004) occurs at much larger distances. If the $\gamma$-ray-emitting jet precesses, the obtained parameters correspond to the average orbital modulation. Still, the observed large modulation amplitude, of $\sim 70–80$ per cent, indicates that the precession does not lead to its substantial reduction.

### 4.3 Modulation of radio emission

We expect to find some orbital modulation of the radio emission in Cyg X-3 caused by free-free absorption in the stellar wind. It is seen, e.g. in the BH HMXB Cyg X-1, where it is strong, with the total amplitude of $\simeq 30 \pm 1$ per cent at 15 GHz (Zdziarski 2012).

We note, however, that the radio observations of Cyg X-3 have been performed with the visibility window repeating each sidereal day, $T_s = 0.99726957$ d. Although the observations were scheduled at times determined by the current collection of other requests for observing, and their priorities, the presence of the visibility window results in a strong peak of the power spectrum around 1 d. In addition, harmonics appear, including the fifth one, which is very close to the orbital period. This has apparently prevented any detection of an orbital modulation of the radio emission in spite of many years of observations available. Indeed, we also do not find a significant peak at the orbital period in the power spectrum of the barycentre-corrected light curve. This is shown in Fig. 10(a), where we see a strong broad peak around 1 d, and the peak around the orbital period is seen at a much lower power. The overall maximum power is at 85.95 d.

Thus, in order to see the orbital modulation in the power spectrum, we follow a technique used in Zdziarski et al. (2012a) for calculating folded light curves. In it, we normalize each flux density to its running linear average, see equation (4) in Zdziarski et al. (2012a), determined in the present case by averaging the flux using the observations within $\pm \delta/2 = 0.1–0.2$ d of its time (i.e. within 1–2 orbital periods) and requiring at least five observations with the positive fluxes in each average. This reduces the number of usable observations by only 4 per cent, and the average number of observational points used for a renormalized flux is $26 \pm 13$ for $\delta = 0.2$ d and $37 \pm 24$ for $\delta = 0.4$ d. We note that this technique corresponds to imposing a high-pass filter in the frequency domain, i.e. it strongly reduces the variability on time-scales longer than the orbital period (very significant in Cyg X-3), thus allowing us to detect the orbital modulation. We also convert the light curve to one corresponding to a constant orbital period in the same way as applied to the $\gamma$-ray light curve, see Section 4.2 above. The power spectrum of the resulting light curve is shown in Fig. 10(b). We see that now the strongest peak is around 0.2 d. We show a zoom of the periodogram to the vicinity of $P_0$ in Fig. 10(c), in which we see that we clearly detect the orbital period of $P_0 = 0.19968476(3)$ of equation (2), with $P_{LS} \simeq 130$.

On the other hand, we also see that the strongest peak of the periodogram is at $P_{LS} \simeq 390$ at a period of $\simeq 0.19944848$ d, clearly different from $P_0$. We have searched for the origin of that peak and found that it corresponds to the fifth harmonic of the sidereal day. In order to clearly see it, we have considered the 15-GHz light curve
<figure><img src="image_12.png" /><figcaption>**Figure 10.** (a) The Lomb–Scargle power spectrum (normalized as in Press et al. 1992) of the barycentre-corrected 15-GHz light curve from the Ryle and AMI telescopes obtained for $\ln F_\nu$. (b) The power spectrum after renormalizing $F_\nu$ to their running averages over ±0.1 d and transforming the time axis to that corresponding to the constant orbital period. (c) The same as in (b) but zoomed to the region containing both the orbital period, $P_0$ (marked by the dashed line), and the strongest peak.</figcaption></figure>

without the barycentric and $\dot{P}$ corrections. However, we still imposed our high-pass filter with $\pm\delta/2 = 0.1$, in order to see variability on time-scales comparable to the orbital period. We have performed the timing analysis on this light curve using both the periodogram and the period-dispersion minimization (PDM; Stellingwerf 1978)

<figure><img src="image_13.png" /><figcaption>**Figure 11.** (a) The results of the PDM analysis of the 15-GHz light curve without both the barycentric correction and the correction for the period increase but after renormalizing $F_\nu$ to their running averages over ±0.1 d, shown for the range of trial periods of 0.1–1.25 d. (b) The same as in (a) but zoomed to the region containing both the orbital period, $P_0$, and the 1/5th of the sidereal day, $T_s$, marked by the dashed and dotted line, respectively. The seen displacement of the former minimum from the exact value of $P_0$ is due to the correction for $\dot{P} > 0$ not being included in this calculation.</figcaption></figure>

method. We show here the results only for the PDM analysis, with those from the periodogram being completely consistent with the PDM ones. Among others, we have found distinct minima of the PDM statistic, $\theta_{\text{PDM}}$, at $T_s$, $T_s/2$, $T_s/3$, $T_s/4$, and $T_s/5$, as shown in Fig. 11(a). A zoom to the region of $T_s/5$ and $P_0$ is shown in Fig. 11(b). We see that the strongest peak seen in Fig. 10(c) corresponds (after removing the time corrections) exactly to $T_s/5$. We have also checked that the period observed in HE γ-rays is equal to $P_0$ and that no additional peaks appear in its vicinity, as shown in Fig. 12, which is consistent with the γ-ray observations being performed from space, thus not affected by the daily visibility window.

Therefore, we hereafter consider only the orbital period. We calculate the orbital modulation of the 15-GHz flux by using the light curve renormalized to the running average and take account of the $\dot{P}$, as described above. We calculate the average flux and its standard deviation within a given phase bin. We first present our results for the entire (renormalized) radio light curve, which corresponds to the modulation averaged over all flux and spectral states of Cyg X-3, see Fig. 13(a). We see a distinct modulation with a depth of ≈4 per cent.

We then separate the data into subsets based on the radio flux and the X-ray spectral state. Based on fig. 2(b) in ZSP16 and Figs 2(b–c) above, we see that the hard state in Cyg X-3 is characterized by the
<figure><img src="image_14.png" /><figcaption>**Figure 12.** The Lomb–Scargle periodogram for Cyg X-3 in the flaring state in the 0.1–100-GeV range (as in Fig. 6, which takes into account the correction for $\dot{P}$) zoomed to the region containing the orbital period (shown by the dashed line) and the strongest peak in Fig. 10(c) (shown by the dotted line).</figcaption></figure>

radio flux density changing in the range of 30 mJy $\lesssim F_v \lesssim 300$ mJy and the 3–5 keV flux of $\lesssim 0.5$ keV cm$^{-2}$ s$^{-1}$. Thus, in order to select the hard state in our data, we use these two criteria for the radio data with available ASM coverage within a day only. Then, the combined soft/intermediate state has the 3–5 keV flux of $\gtrsim 0.5$ keV cm$^{-2}$ s$^{-1}$ and the values of the radio flux density anywhere between undetectable flux and 20 Jy, i.e. also including the 30–300 mJy range. We find that there are virtually no measurements with the 3–5 keV flux being $<0.5$ keV cm$^{-2}$ s$^{-1}$ outside the 30–300 mJy range, so we do not need to impose any condition on that flux. Thus, in order to select the soft/intermediate state in our data, we use the radio data with available ASM coverage within a day with the 3–5 keV flux $>0.5$ keV cm$^{-2}$ s$^{-1}$ for the radio data within 30–300 mJy and all the data with $F_v < 30$ mJy and $F_v > 300$ mJy.

We show the results for the soft/intermediate state in Figs 13(b–d). We see a significant dependence on $F_v$. The modulation amplitude is highest for the lowest radio fluxes, with the amplitude of $\simeq 10$ per cent, and the modulation extrema are at the phases similar to those for the entire data. Then the intermediate and high fluxes show the amplitudes similar to those for the entire data, $\simeq 4$ per cent. There is also a clear shift of the phase of the modulation minima at $\phi/2\pi \simeq 0.6$, $\simeq 1.1–1.3$, and $1.3–1.4$ for the lowest, intermediate, and highest radio-flux range, respectively. Next, we consider the hard state. We also see a significant dependence on $F_v$, see Figs 13(e–f). The modulation amplitude is higher for the lower range of the flux, 30–100 mJy, with the amplitude of $\simeq 6$ per cent, than for the upper range, 100–300 mJy, with the amplitude of $\simeq 2.5$ per cent. We also see a shift of the phase of the modulation minima, at $\phi/2\pi \simeq 0.6$, $\simeq 0.65–0.85$, for the lower and higher radio-flux range, respectively.

The amplitude decrease with the increasing flux within the hard state can be explained by an increase of the height along the jet to where the partially self-absorbed emission becomes optically thin. In the soft/intermediate state, the synchrotron emission is mostly optically thin, but still a larger fraction of the radio flux appears to be emitted close to the compact object for low radio fluxes. Also, the phase of the modulation minimum increasing with the radio flux can be explained by the distance of the location of the bulk of radio emission increasing with the increasing flux and the jet being

<figure><img src="image_15.png" /><figcaption>**Figure 13.** (a) Radio flux folded (using $\ln F_v$) on the orbital period for the entire 15-GHz data renormalized to the local running average (76129). Hereafter, the numbers in parentheses give the respective number of the measurements used. Below, we show the dependences of the orbital modulation on the radio flux range in the soft/intermediate state (b–d) and the hard (e–f) state. We show the results for the data for (b) the lowest soft state ($F_v < 30$ mJy; 5724); (c) the medium soft state ($30$ mJy $< F_v \le 300$ mJy, $F(3–5$ keV) $>0.5$ keV cm$^{-2}$ s$^{-1}$; 14607); (d) the highest soft state ($F_v > 300$ mJy; 8487); (e) the lower hard state ($30$ mJy $< F_v \le 100$ mJy, $F(3–5$ keV) $< 0.5$ keV cm$^{-2}$ s$^{-1}$; 8278); and (f) the upper hard state ($100$ mJy $< F_v \le 300$ mJy, $F(3–5$ keV) $< 0.5$ keV cm$^{-2}$ s$^{-1}$; 8606). The best fit with the wind-absorption model is shown by the solid line in panel (d).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 14.** The mutual dependences of the acceptable parameters (within 68 per cent confidence region), $\phi_j$ versus $\theta_j$, and $H/d$ versus $\beta$, obtained with the isotropic wind model applied to the observed 15-GHz orbital modulation. The colour changes from white to dark blue $\propto \sqrt{n}$, as in Fig. 9.</figcaption></figure>

inclined with respect to the orbital axis. Similar effects are seen in Cyg X-1, see Zdziarski (2012).

The radio orbital modulation and its dependence on the flux can be fitted by the same method and geometry as for $\gamma$-rays (Section 4.2) except that now free-free absorption on the wind is included as the modulation process. We assume a constant wind velocity and temperature, $v_w = 1600$ km s$^{-1}$, $T_w = 10^5$ K, respectively, the mass-loss rate by the donor of $\dot{M}_w = 7.5 \times 10^{-6} M_{\odot} \text{yr}^{-1}$ (Zdziarski et al. 2013 and references therein), the He composition ($X = 0$) and $M_*$, $M_c$, and $i$ as in Section 4.2. Since the radio emission zone is far from the system, we take into account the effect of the non-zero orbital velocity and finite jet velocity, which can cause a significant change in jet orientation at large distances and a phase delay in the radio modulation. These effects are negligible when fitting the $\gamma$-ray orbital modulation.

We fit the orbital modulation averaged over parts of the entire light curve corresponding to the brightest part of the soft state, $F > 0.3$ Jy (Fig. 13d), since the strong $\gamma$-ray emission is observed predominantly in that state and we wish to compare to the jet parameters derived from the $\gamma$-ray modulation. We find that the acceptable regions are wide and include values for the jet inclination, $\theta_j$, and orientation, $\phi_j$, that are compatible with those found for $\gamma$-rays, see Fig. 14. Only the location of the 15-GHz source is robustly constrained (by the amplitude of the modulation) to around $\sim 200d$.

The best-fitting solution is plotted in Fig. 13(d). The region of dominant 15-GHz emission is located at the distance along the jet of $H \simeq 4.4 \times 10^{13}$ cm, for a jet velocity of $\beta \simeq 0.45$, with $\theta_j \simeq 64^{\circ}$, $\phi_j \simeq 140^{\circ}$, and a total $\chi^2 \simeq 1.9$ (10 orbital flux measurements, five fitted parameters). For the best fit, the relative contribution of the counterjet is quite large in the range of 0.33–0.48. A problem with this solution is that it is heavily attenuated, with only a fraction $\simeq 5 \times 10^{-5}$ of the intrinsic flux making it to the observer. We note that the average optical depth, $\langle \tau \rangle$, from the source to infinity is larger by a factor of the order of $\sim H/d$ than the difference in the optical depths between their maximum and minimum values, $\Delta\tau$ (roughly equal to the modulation amplitude), see equations (21–22) in Zdziarski (2012) derived for a perpendicular jet. This explains the large attenuation of this solution.

If we impose the same jet velocity and orientation for both the $\gamma$-ray and radio modulation models, we obtain $\beta \simeq 0.55$, $\theta_j = 30^{\circ}$, and $\phi_j = 43^{\circ}$, with a total $\chi^2 \simeq 9.5$ (with contributions of 5.5 and 4.0 for the $\gamma$-ray and radio part, respectively). The locations of the sources are $H_\gamma \simeq 1.1 \times 10^{12}$ cm and $H_r \simeq 4.7 \times 10^{13}$ cm. The average attenuation of this solution is more moderate, $5.4 \times 10^{-2}$.

Finally, we mention that in our investigations we also considered the hypothesis that the strongest peak in Fig. 10(c) is due to a beat with a precession of the jet. The difference of the period of the strongest peak and the orbital period of $\simeq 20$ s corresponds to a precession period of $\simeq 170$ d (which, interestingly, is about twice the period of the strongest peak in the periodogram of Fig. 10a at $\simeq 86$ d). We have then searched for a dependence of the orbital modulation on the precession phase, and, surprisingly, we have found a rather regular dependence on it of both the precession amplitude and the orbital phases of the modulation extrema. Still, the exact coincidence of the found period with the fifth harmonic of the sidereal day convinced us of its origin as an artefact of the visibility window of the radio telescopes. The regular behaviour we found could thus be spurious.

# 5 CROSS-CORRELATIONS

## 5.1 The method

We calculate the Pearson's correlation coefficient between a discrete light curve, $x_i$, $i = 1, \dots, I$, and another one, $y_j$, $j = 1, \dots, J$, shifted in time by $\Delta t$,

$$r(\Delta t) = \frac{\sum [x_i - \bar{x}(\Delta t)][y_i - \bar{y}(\Delta t)]/K}{\sqrt{\sum_i [x_i - \bar{x}(\Delta t)]^2/I'} \sqrt{\sum_j [y_j - \bar{y}(\Delta t)]^2/J'}}, \quad (3)$$

where the summation in the numerator is over all pairs, $(i, j)$, satisfying

$$\Delta t - \delta/2 \le t(y_j) - t(x_i) < \Delta t + \delta/2, \quad (4)$$

$\delta$ is the bin size of the time shift, $\Delta t$, $K$ is the number of such pairs, and the values of $\bar{x}$ and $\bar{y}$ and the sums in the denominator and in $\bar{x}$, $\bar{y}$, are over only $I'$, $J'$, values of the $x_i$ or $y_j$ satisfying equation (4), respectively. The standard deviation of $r$ is calculated using equation (5) of Edelson & Krolik (1988).

This method differs slightly from that of Edelson & Krolik (1988), who considered the values of the standard deviations and the averages for each light curve based on all of their respective points, while here we include only those entering a given bin, as proposed by Lehár et al. (1992). Using global averages and standard deviations can lead to substantial inaccuracies if either a light curve has long-term trends and the cross-correlation is carried over
a section of it or it is strongly varying. In particular, we found that the value of the auto-correlation at zero lag is substantially greater than unity in a number of cases considered here. In fact, the value of *r* is within the range of [-1, 1] only if *I* = *J* = *K*. Imposing that requires multiple counting (for each occurrence of condition 4) in the mean and standard deviation in cases in which a given *x*ᵢ satisfies condition (4) for more than one *y*ᵢ (or vice versa). Similarly to the case of using global averages and standard deviations, the auto-correlation can exceed unity and the cross-correlation can be not correctly normalized if the multiple counting for the mean and standard deviation is not allowed.

Then, as an option, we average each light curve in the pair within its bin of the size δ before calculating *r*. This alleviates the above problem of the normalization of *r*, since then a given *x*ᵢ satisfies condition (4) for (typically) only one value of *y*ᵢ. We have found this to be especially important in the case of correlating the 15-GHz light curve from the AMI with the LAT γ-rays, in which case the LAT light curve has one point per day while the AMI one has typically several. We find then the correlation coefficient with a very noisy dependence on Δ*t*, caused by variations of both the average values and the standard deviations. However, for each shown correlation, we have tested that using different options leads to similar overall shapes of the correlations. Also, we use logarithms of the fluxes, since the flux distributions in Cyg X-3 are much closer to lognormal than to normal (ZSP16), and the calculation of *r* of equation (3) assumes that the distributions of *x*ᵢ and *y*ᵢ are normal. Lognormal flux distributions have been found in other accreting systems (Uttley, McHardy & Vaughan 2005).

## 5.2 Cross and auto-correlations in Cyg X-3

We first cross-correlate the γ-ray and radio light curves. We show the results for the LAT detections with the fractional error ≤0.5 within a day (yielding 174 d) in Fig. 15(a) (hereafter the dashed curves show the uncertainty ranges estimated as above). The highest γ-ray fluxes are found in the soft and intermediate states, though we also find a relatively large number of detections in the hard state. The cross-correlation peaks at zero lag, implying that the lag averaged over all frequencies of the variability is <1 d. This is a much more accurate result than the early one of FLC09, who obtained a peak lag of 5 ± 7 d. However, the cross-correlation shows a significant asymmetry, indicating that some radio photons still lag the γ-ray emission by ≲10 d. When the required maximum fractional error is increased, the cross-correlation still peaks at zero lag but becomes weaker, with a lower value of *r*. We have also considered the case with the radio data split into parts with $F_ν > 0.3$ and <0.3 Jy. For both ranges, the cross-correlations peak at zero lag. Fig. 16(a) shows the relationship between the two daily-averaged energy fluxes on the same MJDs, where we see a positive correlation of the γ-ray flux with the radio one, $F_γ \propto F_R^{1/3}$.

The γ-ray emission (all detections with the fractional error <1; 486 d) has a narrow auto-correlation with the width of <1 d, see Fig. 15(b). We also see a weak auto-correlation tail, dropping to null at ≳10 d. On the other hand, the radio emission has a relatively wide auto-correlation, with the half-width at *r* = 0.5 of 4.5 d, see Fig. 15(c). We then see that the radio emission is anti-correlated with itself for Δ*t* ~ 30 d. In order to investigate the origin of it, we have split the data set into two parts, above and below 0.3 Jy. We find that the auto-correlation for the low fluxes is similar to the one for all the data, and the ~30 d time-scale appears to be related to the typical duration of a single occurrence of the hard state. On the other hand, the auto-correlation for the high fluxes is

<figure><img src="image_17.png" /><figcaption>**Figure 15.** (a) The cross-correlation (solid curve) between the 1-d LAT detections with the fractional error ≤0.5 and the 15-GHz radio emission from the AMI. Hereafter, Δ*t* > 0 corresponds to the signal in the second photon energy range given in the plot label lagging behind the signal in the first range (in the present case, 15-GHz flux lagging the γ-ray one), and dashed curves give the estimated uncertainty range of a correlation. Δ*t* > 0 corresponds to the radio emission lagging the γ-rays. We do not find a measurable lag of the radio emission. (b) The auto-correlation of the γ-ray detections (with the fractional error <1 within a day). (c) The auto-correlation of the 15-GHz emission during the same epoch as that of the Fermi observations.</figcaption></figure>

narrower and it becomes negative already at Δ*t* ≃ 8 d and reaches the global negative minimum at 20 d. This appears to be related to the radio flares both preceded and followed by states with weak radio emission (e.g. Szostek et al. 2008).

We then correlate the X-ray and γ-ray emission. Fig. 16(b) shows the correlation with the soft X-rays on the same MJDs. We see that almost all γ-ray detections in the hard state have relatively low fluxes, below the boundary of the flaring state. Then, there is a positive correlation in the soft state. The cross-correlation with soft X-rays is shown in Figs 17(a–b), where the positive correlation at zero lag continues up to a few tens of days, indicating that the γ-ray emission continues after an occurrence of a peak in soft X-rays.
<figure><img src="image_18.png" /><figcaption>**Figure 16.** The relationships between the daily 0.1–100 GeV energy fluxes with the fractional error <0.5 and the average of (a) the 15-GHz radio flux, (b) the ASM 3–5 keV flux, and (c) the BAT 15–50 keV flux, on the corresponding MJDs, respectively. The horizontal dotted line gives the boundary of the flaring γ-ray state, and the vertical dashed lines give the boundaries between the hard and soft/intermediate states. In (a), the hard state occurs only between the two dashed lines, but this range also corresponds to the soft/intermediate state. In (b) and (c), the soft state occurs to the right and left of the dashed line, respectively.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Left-hand panels: The cross-correlation between (a) the 3–5 keV ASM rate and (b) 2–4 keV MAXI rate versus the LAT detections (positive lags correspond to γ-rays lagging the X-rays). Middle panels: The cross-correlation between the 3–5 keV ASM rate in (c) the hard state ($F_X < 0.5$ keV cm$^{-2}$ s$^{-1}$, 30 < $F_ν$ < 300 mJy) and in (d) the soft/intermediate state ($F_X > 0.5$ keV cm$^{-2}$ s$^{-1}$). Right-hand panels: The cross-correlation between the 15–50 keV BAT rate in (e) the hard state ($F_X > 1.1$ keV cm$^{-2}$ s$^{-1}$, 30 < $F_ν$ < 300 mJy) and in (f) the soft/intermediate state ($F_X < 1.0$ keV cm$^{-2}$ s$^{-1}$). Positive lags correspond to radio lagging the X-rays. For all data sets, only data with the fractional error $\le 0.5$ were considered.</figcaption></figure>

This corresponds to occurrences of γ-ray flares during transitions from the softest (hypersoft) X-ray states to harder ones, confirming Koljonen et al. (2010). However, some γ-ray flares also take place during the opposite transitions, which may correspond to weaker peaks at negative lag in Figs 17(a–b),

The correlation with hard X-rays is more complicated, see Fig. 16(c). While most of γ-ray detections with high fluxes occur for low hard X-ray flux, i.e. in the soft state, which corresponds to an overall anti-correlation, the detections within the soft state show a weak positive correlation. We also see a large number of detections in the hard state below the boundary of the flaring state. This results in a strong anti-correlation between the γ-ray emission and hard X-rays, with the minimum of the cross-correlation at a lag of γ-rays of ~5 d (not shown here). However, that lag, given the measurement errors, may be not statistically significant.

We next consider cross-correlations between X-rays and radio in different states. We show the cross-correlations between the soft X-rays (3–5 keV) and radio in the hard and soft/intermediate states in Figs 17(c) and (d), respectively. We see that the hard-state relationship gives a strong positive correlation peaking at zero lag (for 1-d bins). The cross-correlation is relatively wide, with a half-width of ~15 d, and relatively symmetric, showing that some radio photons in the hard state lead and some lag the soft X-rays. We also see some asymmetry at $|\Delta t| > 20$ d, indicating that the lag dominates at long $\Delta t$.

On the other hand, while the soft/intermediate state shows almost no correlation at zero lag, it shows a strong positive peak at the radio band delayed by $\approx 45–48$ d, see Fig. 17(d). To investigate it further, we have calculated the relationship between the soft X-ray and radio fluxes in the soft/intermediate state at 0- and 46-d lags. Those
<figure><img src="image_20.png" /><figcaption>**Figure 18.** (a) The auto-correlation for the 3–5-keV ASM (upper blue curve) and the 15–50-keV BAT (lower red curve) rates. (b) The cross-correlation between the 3–5-keV and 15–50-keV rates in all states.</figcaption></figure>

plots, not shown here, confirm the lack of a correlation at zero lag changing into a weak positive correlation at the 46-d lag. We note that the 45–48-d lag appears to correspond to the average time spent in the ultra/hypersoft X-ray state that directly precedes major radio flares. For instance, in 2011, the radio emission was quenched for a month while the source was in the ultra/hypersoft state, ending with a major 10-Jy radio flare (Corbel et al. 2012). The interpretation of the physical nature of this lag appears, however, unclear. It may correspond to a time-scale linked to magnetic field rearrangement in the disc so that a jet can be launched. It may also correspond to the propagation time-scale from the centre to the dominant large-scale jet component of the radio emission (at tens of mas; Tudose et al. 2010). This interpretation, however, implies a rather low speed of the jet, ~0.1c. The radio emission is also seen on arcsec scales, as found by Martí et al. (2001), which may contribute to that long lag as well. We have also calculated the soft X-rays versus 15-GHz cross-correlation without separating into the states to be able to directly compare it to the 3–5 keV versus γ-ray cross-correlation. It shows the shape relatively similar to that of Fig. 17(b), with a positive *r* at zero lag and a peak below 50 d.

Hard X-rays, 15–50 keV, are anti-correlated with the radio band in the hard state (i.e. for large X-ray fluxes), see Fig. 17(e). There appears to be a 1-d lag of the radio emission here, but it is not statistically significant. On the other hand, the hard X-rays are positively correlated with the radio at zero lag in the soft/intermediate state, see Fig. 17(f). However, they also show an anti-correlation with radio peaking at a lag of ~30 d; see Fig. 17(f). This behaviour is likely to be related to the 45–48-d lag at soft X-rays, Fig. 17(d).

We then show the auto-correlation functions for the 3–5 and 15–50 keV energy ranges, and their cross-correlations in Figs 18(a–b), respectively. Interestingly, the 15–50-keV auto-correlation is substantially narrower, with the half-width of ≈20 d, than the 3–5 keV one, with the half-width of ≈40 d. This difference may

correspond to the hard X-ray emission region being closer to the compact object and thus smaller than that for soft X-rays. On the other hand, it may also correspond to the 15–50-keV emission leading the 3–5 keV one. Interestingly, the widths of the X-ray auto-correlations are of the same order as the radio versus X-ray lags in the soft state, Figs 17(d–f).

We also note that Tudose et al. (2010) argued that the radio/X-ray correlation observed at zero lag (e.g. Szostek et al. 2008) is not theoretically expected in bright radio states, where the bulk of the radio emission is in the jet rather than in the core. However, this does not seem to present a problem, given the projected jet distance from the core of ~1–2 light days and the widths of the auto-correlations of both radio and X-rays are ≫2 d.

## 6 CONCLUSIONS

We have obtained the following main results.

Based on 9 yr of the *Fermi* data, we have searched for occurrences of significant HE γ-ray emission from Cyg X-3. We have found a large number of days with significant LAT detections, and among them, 49 d with the 0.1–100-GeV energy fluxes above the flaring level, which we defined at $4.94 \times 10^{-10} \text{ erg cm}^{-2} \text{ s}^{-1}$ based on the level $3\sigma$ above the average of non-flaring detections. Out of them, 43 d are in the soft spectral state.

We have calculated the average γ-ray spectrum during strong flares (with positive detections in the 0.08–15-GeV energy range) and the spectrum averaged over all the occurrences of the soft spectral state (detected in the 0.5–15-GeV range). On the other hand, we have found only upper limits in the hard and intermediate states. The flaring-state spectrum is well modelled by Compton scattering of the blackbody photons from the donor by jet relativistic electrons with a power-law distribution with the spectral index of ≈3.5 and the low-energy cutoff at the Lorentz factor of ~10³.

From the LAT data, we have also obtained the profile of the orbital modulation of γ-rays in the flaring state, which is significantly more accurate than the previous one of FLC09. The amplitude of the modulation is large, by a factor of ≈5, and it is well modelled by Compton scattering of stellar blackbody, which agrees with the modelling of the spectrum. The modulation model implies the location of the γ-ray source at the distance along the jet similar to the separation between the binary components, a mildly relativistic jet velocity, and it requires that the jet is inclined at an angle ≳25° with respect to the binary axis.

We have then studied the 22 yr of 15-GHz radio observations of Cyg X-3 by the Ryle and AMI telescopes. We have discovered pronounced modulation of the radio emission at the orbital period. The amplitude of the modulation depends on both the spectral state and the flux level. It changes from ≈2.5 to ≈10 per cent, and it is ≈4 per cent when averaging over all the data. We model the observed modulation as free–free absorption in the stellar wind of the jet radio emission. We find the 15-GHz source to be located at a distance of ~10² times the binary separation.

Finally, we have studied cross-correlations between the HE γ-ray and radio light curves, as well as between either of them and light curves from the ASM, MAXI, and BAT X-ray monitors. We have found the correlation coefficient between the HE γ-ray and 15-GHz light curves peaks at zero lag. However, its asymmetry indicates that some radio photons lag the γ-rays by <10 d. Then, the γ-rays lag soft X-rays by some tens of days but without showing a clear peak of the correlation coefficient. This is consistent with occurrence of γ-ray flares mostly during the soft-to-hard transitions. Also, we have not found measurable lags between the X-ray and radio emission
in the hard spectral state. On the other hand, we have found the lag peaking at 45–48 d of the 15-GHz emission with respect to 3–5-keV soft X-rays in the soft state. This lag is similar to the time spent in the X-ray ultra/hypersoft and radio-quenched state before a major flare occurs.

## ACKNOWLEDGEMENTS

We thank the referee and members of the *Fermi* LAT Collaboration for valuable comments. This research has been supported in part by the Polish National Science Centre grants 2013/10/M/ST9/00729, 2015/18/A/ST9/00746, and 2016/21/P/ST9/04025, and the Carl-Zeiss Stiftung through the grant 'Hochsensitive Nachweistechnik zur Erforschung des unsichtbaren Universums' to the Kepler Zentrum für Astro- und Teilchenphysik at the University of Tübingen, and by the state of Baden-Württemberg through bwHPC. The authors thank the SFI/HEA Irish Centre for High-End Computing (ICHEC) for the provision of computational facilities and support. The *Fermi* LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work has been performed in part under DOE Contract DE-AC02-76SF00515. The AMI arrays are supported by the University of Cambridge and by the UK STFC. This research has made use of the results provided by the ASM and BAT teams and MAXI data provided by RIKEN, JAXA, and the MAXI team.

## REFERENCES

Ackermann M. et al., 2012, ApJ, 750, 3

Atwood W. B. et al., 2009, ApJ, 697, 1071

Axelsson M., Larsson S., Hjalmarsdotter L., 2009, MNRAS, 394, 1544

Barthelmy S. D. et al., 2005, Space Sci. Rev., 120, 143

Belczyński K., Bulik T., Mandel I., Sathyaprabakar B. S., Zdziarski A. A., Mikołajewska J., 2013, ApJ, 764, 96

Bhargava Y. et al., 2017, ApJ, 849, 141

Bodaghee A., Tomsick J. A., Pottschmidt K., Rodriguez J., Wilms J., Pooley G. G., 2013, ApJ, 775, 98

Bradt H. V., Rothschild R. E., Swank J. H., 1993, A&AS, 97, 355

Bulgarelli A. et al., 2012, A&A, 538, A63

Burke M. J., Gilfanov M., Sunyaev R., 2017, MNRAS, 466, 194

Cheung C. C., Loh A., 2016, Astron. Telegram, 9502

Corbel S., Fender R. P., Tzioumis A. K., Nowak M., McIntyre V., Durouchoux P., Sood R., 2000, A&A, 359, 251

Corbel S. et al., 2012, MNRAS, 421, 2947

Corbel S., Coriat M., Brocksopp C., Tzioumis A. K., Fender R. P., Tomsick J. A., Buxton M. M., Bailyn C. D., 2013, MNRAS, 428, 2500

Dubus G., Cerutti B., Henri G., 2010, MNRAS, 404, L55

Edelson R. A., Krolik J. H., 1988, ApJ, 333, 646

Eggleton P. P., 1983, ApJ, 268, 368

Egron E. et al., 2017, MNRAS, 471, 2703

Fender R. P., Hanson M. M., Pooley G. G., 1999, MNRAS, 308, 473

Fender R. P., Belloni T. M., Gallo E., 2004, MNRAS, 355, 1105

Fermi-LAT Collaboration, 2009, Science, 326, 1512 (FLC09)

Fermi-LAT Collaboration, 2010, ApJS, 188, 405

Fermi-LAT Collaboration, 2015, ApJS, 218, 23 (3FGL)

Giacconi R., Gorenstein P., Gursky H., Waters J. R., 1967, ApJ, 148, L119

Hanson M. M., Still M. D., Fender R. P., 2000, ApJ, 541, 308

Hjalmarsdotter L., Zdziarski A. A., Larsson S., Beckmann V., McCollough M., Hannikainen D. C., Vilhu O., 2008, MNRAS, 384, 278

Hjalmarsdotter L., Zdziarski A. A., Szostek A., Hannikainen D. C., 2009, MNRAS, 392, 251

Koljonen K. I. I., Maccarone T. J., 2017, MNRAS, 472, 2181

Koljonen K. I. I., Hannikainen D. C., McCollough M. L., Pooley G. G., Trushkin S. A., 2010, MNRAS, 406, 307

Koljonen K. I. I., Maccarone T., McCollough M. L., Gurwell M., Trushkin S. A., Pooley G. G., Piano G., Tavani M., 2018, A&A, 612, A27

Krimm H. A. et al., 2013, ApJS, 209, 14

Lehár J., Hewitt J. N., Roberts D. H., Burke B. F., 1992, ApJ, 384, 453

Levine A. M., Bradt H., Cui W., Jernigan J. G., Morgan E. H., Remillard R., Shirey R. E., Smith D. A., 1996, ApJ, 469, L33

Ling Z., Zhang S., Tang S., 2009, ApJ, 695, 1111

Loh A., Corbel S., 2017a, Astron. Telegram, 10109

Loh A., Corbel S., 2017b, Astron. Telegram, 10243

Loh A. et al., 2016a, MNRAS, 462, L111

Loh A., Corbel S., Dubus G., Corbet R., 2016b, Astron. Telegram, 8591

Lyne A. G., Stappers B. W., Keith M. J., Ray P. S., Kerr M., Camilo F., Johnson T. J., 2015, MNRAS, 451, 581

Markwardt C. B., Tueller J., Skinner G. K., Gehrels N., Barthelmy S. D., Mushotzky R. F., 2005, ApJ, 633, L77

Martí J., Paredes J. M., Peracaula M., 2001, A&A, 375, 476

Matsuoka M. et al., 2009, PASI, 61, 999

Mattox J. R. et al., 1996, ApJ, 461, 396

McCollough M. L. et al., 1999, ApJ, 517, 951

McCollough M. L., Corrales L., Dunham M. M., 2016, ApJ, 830, L36

Miller-Jones J. C. A., Blundell K. M., Rupen M. P., Mioduszewski A. J., Duffy P., Beasley A. J., 2004, ApJ, 600, 368

Mioduszewski A. J., Rupen M. P., Hjellming R. M., Pooley G. G., Waltman E. B., 2001, ApJ, 553, 766

Munar-Adrover P., Paredes J. M., Romero G. E., 2011, A&A, 530, A72

Piano G. et al., 2012, A&A, 545, A110

Piano G. et al., 2016, Astron. Telegram, 9429

Piano G. et al., 2017a, Astron. Telegram, 10138

Piano G. et al., 2017b, The Astron. Telegram, 10179

Pooley G. G., Fender R. P., 1997, MNRAS, 292, 925

Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical Recipes in FORTRAN. The Art of Scientific Computing. Cambridge University Press

Stellingwerf R. F., 1978, ApJ, 224, 953

Szostek A., Zdziarski A. A., 2008, MNRAS, 386, 593

Szostek A., Zdziarski A. A., McCollough M. L., 2008, MNRAS, 388, 1001

Takata J., Tam P. H. T., Ng C. W., Li K. L., Kong A. K. H., Hui C. Y., Cheng K. S., 2017, ApJ, 836, 241

Tavani M. et al., 2009a, A&A, 502, 995

Tavani M. et al., 2009b, Nature, 462, 620

Tudose V. et al., 2007, MNRAS, 375, L11

Tudose V. et al., 2010, MNRAS, 401, 890

Uttley P., McHardy I. M., Vaughan S., 2005, MNRAS, 359, 345

van der Klis M., Bonnet-Bidaud J. M., 1989, A&A, 214, 203

van Kerkwijk M. H., 1993, A&A, 276, L9

van Kerkwijk M. H. et al., 1992, Nature, 355, 703

van Kerkwijk M. H., Geballe T. R., King D. L., van der Klis M., van Paradijs J., 1996, A&A, 314, 521

Vilhu O., Hakala P., Hannikainen D. C., McCollough M., Koljonen K., 2009, A&A, 501, 679

Vladimirov A. E. et al., 2011, Comput. Phys. Commun., 182, 1156
Williams P. K. G. et al., 2011, ApJ, 733, L20

Yoon D., Zdziarski A. A., Heinz S., 2016, MNRAS, 456, 3638

Zanin R., Fernández-Barral A., de Oña-Wilhelmi E., Aharonian F., Blanch
O., Bosch-Ramon V., Galindo D., 2016, A&A, 596, A55

Zdziarski A. A., 2012, MNRAS, 422, 1750

Zdziarski A. A., Misra R., Gierliński M., 2010, MNRAS, 402, 767

Zdziarski A. A., Maitra C., Frankowski A., Skinner G. K., Misra R., 2012a,
MNRAS, 426, 1031

Zdziarski A. A., Sikora M., Dubus G., Yuan F., Cerutti B., Ogorzałek A., 2012b, MNRAS, 421, 2956

Zdziarski A. A., Mikołajewska J., Belczyński K., 2013, MNRAS, 429, L104

Zdziarski A. A. et al., 2016a, MNRAS, 455, 1451

Zdziarski A. A., Segreto A., Pooley G. G., 2016b, MNRAS, 456, 775 (Er-atum: 2016, MNRAS, 457, 3887) (ZSP16)

Zdziarski A. A., Malyshev D., Chernyakova M., Pooley G. G., 2017,
MNRAS, 471, 3657

This paper has been typeset from a TeX/LA TeX file prepared by the author.