© F. Cangemi et al. 2021

1. Introduction

Black-hole X-ray binaries (BHXBs) are systems with variable X-ray emission and can transit through various spectral states during outbursts. Because of the discovery of different types of radio relativistic jets in these systems (e.g., Mirabel et al. 1992; Fender et al. 1999), they are commonly dubbed “microquasars”. The two main (canonical) states are the low-hard state (LHS) and the high-soft state (HSS; see, e.g., Remillard & McClintock 2006; Belloni 2010, for a precise definition of the states). In the LHS, the X-ray spectrum between ∼1 and ∼100 keV is dominated by a hard power law with a photon index of Γ ≲ 1.5 which cuts off at a few hundred keV. Compact radio jets are only observed in the LHS (e.g., Fender 2001) and have been resolved in a couple of sources: Cygnus X-1 (Cyg X-1) and GRS 1915+105 (e.g., Stirling et al. 2001; Fuchs et al. 2003). In the HSS, the spectrum is dominated by a blackbody component peaking at ∼1 keV. This component is associated with thermal emission from an optically thick accretion disk. Usually, no radio emission is detected in this state (e.g., Fender et al. 1999; Remillard & McClintock 2006; Belloni 2010), with the possible exception of Cyg X-1 where tenuous radio emission was detected using very long baseline interferometry (Rushton et al. 2012).

In several bright sources, observations made at higher energies reveal the presence of a high-energy tail above a few hundred keV that can sometimes extend to several MeV (Grove et al. 1998). In this paper we seek further insight into the state-resolved high-energy tail emission in order to discuss its physical origin which is still subject to much debate. To do so, we focus on the most promising candidate, Cyg X-1, and exploit the combined capabilities of the Joint European X-Ray Monitor (JEM-X, Lund et al. 2003), the Imager on Board the INTEGRAL Satellite (IBIS; Lebrun et al. 2003), and the SPectrometer onboard INTEGRAL (SPI; Vedrenne et al. 2003) to perform detailed spectral analysis of the broad-band ∼3 keV–2 MeV data by accumulating all existing INTEGRAL data of this source taken over the first 15 years of the satellite’s life.

Although it was used as a prototype for the definition of the (canonical) X-ray spectral states, Cyg X-1 is a peculiar case. Located at a distance of 1.86 ± 0.12 kpc (Reid et al. 2011), it is a persistent high-mass X-ray binary (HMXB) with a 14.8 ± 1.0 M_⊙ black hole orbiting a blue supergiant star with a mass of 19.2 ± 1.9 M_⊙ (Orosz et al. 2011). The orbital plane is seen at an inclination of ∼27° (Orosz et al. 2011), a value that we use throughout this paper. In contrast to the majority of (transient) BHXBs, its individual states persist on much longer timescales (e.g., Grinberg et al. 2013), and even though the parameters ofx its spectral states are similar to those of other sources, Cyg X-1 does not follow the typical “q-track” pattern for transients: for example, it has never returned to quiescence.

A high-energy tail above the thermal Comptonization continuum in Cyg X-1 was detected for the first time in both the soft and hard accretion states by the COMPTEL instrument (McConnell et al. 2000, 2002) on board the Compton Gamma-Ray Observatory (CGRO). The tail in the LHS was later confirmed with the SPI (Cadolle Bel et al. 2006; Jourdain et al. 2012) and IBIS (Laurent et al. 2011; Rodriguez et al. 2015, hereafter L11, R15), but so far has remained undetected by INTEGRAL in the HSS (Jourdain et al. 2014, hereafter J14). This might be both due to the intrinsic faintness of the component and the rather short exposure time of INTEGRAL in the HSS. Observations obtained between 2003 and 2011 contained ∼90% of LHS data (e.g., R15). Thanks to a spectral transition to a long-lasting HSS in 2011 (Grinberg et al. 2011), follow-up studies such as the present one will be increasingly sensitive to the high-energy tail emission of Cyg X-1 in the HSS.

Detection and investigation of the high-energy tails in both states of Cyg X-1 with INTEGRAL are possible thanks the mission’s combination of a large effective area with the capability to measure polarization – a key tool for breaking degeneracies in the quest for the physical origin of the high-energy tail emission. Several studies have confirmed a high degree of polarization for the high-energy tail in the LHS (L11, Jourdain et al. 2012, and R15 for a state-resolved polarimetric analysis), which is compatible with synchrotron emission from a jet (observed in the source’s LHS e.g., Stirling et al. 2001; Zdziarski et al. 2012). However, this extension above the Comptonized spectrum could also originate from a nonthermal electron distribution in the corona (e.g., McConnell et al. 2002; Romero et al. 2014; Jourdain et al. 2014; Chauvin et al. 2018). In particular, Del Santo et al. (2013) separated all Cyg X-1 INTEGRAL data obtained until 2009 into 12 separate sub-states according to their hardness in the IBIS range and their 30–100 keV spectral slope. They studied these 12 stacked spectra within the framework of hybrid models, either nonmagnetized (eqpair) or with magnetization (belm), and were able to obtain upper limits on the corona’s magnetic field.

We are therefore not only interested in assessing the existence of a high-energy tail in the HSS but also in investigating its origin with the extended amount of data available. Our approach, accumulating data until 2017, increases the count statistics by a significant amount over all earlier works and allows us to study the parameter space of hybrid models not only in the LHS but also, for the first time, with INTEGRAL in the HSS up to 2 MeV. This paper is organized as follows. In Sect. 2, we outline the selection criteria for data considered in this work and the extraction of data products. In Sect. 3, we present a set of techniques and models used to describe the data products. In Sect. 4, we discuss the results and in Sect. 5 we provide our conclusions.

2. Observations and data reduction

2.1. State classification and data selection

We consider all Cyg X-1 observations made with INTEGRAL from 2003 June 9 to 2017 November 7 (Modified Julian Day, MJD 52799 to 58054) corresponding to the satellite revolutions (hereafter revs) 54 to 1882. We restrict our analysis to “science windows” (hereafter scws), that is, individual INTEGRAL pointings of typically 1800–3600 s duration, where Cyg X-1 is less than 10° off axis. This results in a total of 7288 scws corresponding to 11.4 Ms of dead time-corrected exposure.

Similar to R15, we separate the whole data set into three states (soft/HSS, intermediate/IS, and hard/LHS) using the model-independent classification of Grinberg et al. (2013). This classification defines LHS, HSS, and IS based on count rates and hardnesses measured by all-sky monitor instruments (RXTE/ASM, MAXI, and Swift/BAT). Figure 1 shows the classification of the individual scws. This figure extends the similar Fig. 1 from Grinberg et al. (2013), to whom we refer for the details of the methodology of the classification. As in R15, we use all-sky monitor data within 6 hours from a given INTEGRAL scw for the classification. Such stringent time constraints on the classification are necessary because the source is known to undergo fast state transitions (Böck et al. 2011; Grinberg et al. 2013). In total, we roughly double the total exposure time in the LHS compared to R15, and the new extended HSS covering roughly MJD 55800 to 57400 (with some excursions into the LHS, Figure 1) allows us to increase the HSS exposure by a factor of three compared to R15. After MJD 57400, INTEGRAL mostly catches Cyg X-1 in the LHS (Fig. 1). To summarize, we obtain 3178, 2262, and 600 scws, in the LHS, HSS, and IS respectively, corresponding to 4.4 Ms, 3.4 Ms, and 0.9 Ms of exposure time. A total of 1248 scws remain unclassified with this method and are not taken into account in our analysis.

Fig. 1. Spectral state evolution of Cyg X-1 between 2003 and 2017 after Grinberg et al. (2013). The HSS is shown in red, the IS in green, and the LHS in blue. Unclassified data are plotted in gray. Upper panel: 30–50 keV INTEGRAL/ISGRI light curve. The vertical dashed line indicates the change of the INTEGRAL data analysis software. Lower three panels: RXTE/ASM, MAXI, and Swift/BAT light curves of Cyg X-1.

2.2. Spectral extraction: INTEGRAL/IBIS/ISGRI

Data from INTEGRAL/IBIS’ upper detector layer ISGRI¹ taken before December 2015 (INTEGRAL revolution 1627, MJD 57387) were reduced with version 10.2 of the Off-line Scientific Analysis (OSA) software. Data obtained after that date were extracted with OSA v11. This split is necessary because (at the time of writing) specific versions of calibration data —specifically for the gain correction of the instrumental energy scale— are required for both software versions². The long-term IBIS light curve (with a resolution of individual scws) is shown in the top panel of Fig. 1. The different gain correction used before and after December 2015 is seen as an overall increase in the 30–50 keV count rate. For comparison, the 30–50 keV Crab count rate is ∼87 counts s⁻¹ for OSA v10.2 and ∼134 counts s⁻¹ for OSA v11.

IBIS/ISGRI data were reduced with the standard procedure³. For each satellite revolution we first produce a mosaic image to identify the brightest sources (here defined as those with a detection significance greater than 7σ in the 30–50 keV mosaic image) of the field. We then extract light curves and spectra on a scw-by-scw basis. For each scw, we construct a sky model including the brightest sources as defined above (usually Cyg X- 3, Cyg X-2, GRO J2058+42, and KS 1947+300, in addition to Cyg X-1). Taking the background into account, the analysis software then reconstructs the sky image and corresponding source count rates by deconvolving the shadowgrams projected on the detector plane. The roughly 250 scws from the time period of the 2015 June outburst of V404 Cyg (INTEGRAL revolutions 1554 through 1558) remain unclassified and were therefore ignored in our analysis.

Source spectra were extracted on a logarithmic grid of 60 spectral channels. In OSA v11.0, the spectral binning and the response matrices are calculated and provided automatically by the extraction software. In OSA v10.2, the binning is provided by the user and based on a pre-rebinned redistribution matrix file (rmf). We construct the latter to be as close as possible to the OSA v11.0 one, resulting in spectral bin deviations of less than 0.25 keV per channel. We then stack the individual spectra according to the version of the extraction software and following their spectral state classification. This results in a total of four stacked spectra: two per spectral state. In all spectra, we add 2% systematic uncertainty to each of the channels⁴. For the spectral fittings, we omit the spectral channels below 25 keV for all the ISGRI spectra obtained with OSA v10.2 and those below 30 keV for the OSA 11.0 ones.

2.3. Spectral extraction: INTEGRAL/SPI

The SPI data analysis relies on the comparison of sky models and a description of the instrumental background to data in the raw format of counts per energy bin, detector, and pointing (scw). See Vedrenne et al. (2003) for details about the instruments. For the SPI extraction we include the source-intrinsic and background variability in the sky model, yielding state-resolved, long-term spectra. In the case of Cyg X-1 our sky models consist of point-sources in the 16° ×16° field of view with time-dependent fluxes. Above ∼100 keV, only Cyg X-1 and Cyg X-3 are significantly detected and taken into account.

First, we collect the good time intervals as defined for the IBIS/ISGRI data analysis divided into soft and hard state. Because SPI is very sensitive to instrument mode changes, (solar) particle events, detector failures, and Van Allen belt passages, we apply additional filter criteria as shown in Table 1. The resulting number of pointings is 2358 and 1771 for the LHS and HSS, respectively, with a total dead-time-corrected observation time of 3.4 Ms and 4.7 Ms, respectively. The resulting exposure map using all good pointings is shown in Fig. 2. This map shows Cyg X-1 in the center of the chosen exposure and which (known) sources are exposed and how they might contribute in the fit.

Fig. 2. SPI exposure of Cyg X-1 in its soft state in units of 106 cm2 s−1. The fully coded field of view of SPI is 16° ×16°. Green points indicate point sources that have to be taken into account in the SPI data analysis. Red points show other sources listed in the SPI catalog of Bouchet et al. (2011) but only seen below ≲100 keV or at specific flare times.

We consider single events (SEs) below 530 keV and pulse-shape discriminated (PSD) events above that energy. This removes strong “electronic noise” features between 1.4 and 1.7 MeV and allows for a smooth spectral extraction⁴ We also extract SE spectra up to 755 keV to establish consistency between the two data types. The derived PSD fluxes are then scaled by a factor of 1/0.85 to account for the change in efficiency ⁴.

For the spectral extraction, we define 21 logarithmic energy bins between 25 keV and 755 keV for the SE range, and 8 logarithmic bins between 530 keV and 2000 keV for the PSD range. The spectra are obtained bin by bin from maximum-likelihood fits accounting for Poisson-distributed data. These fits include the two-component background model as described by Siegert et al. (2019), which is based on the SPI BackGround and Response DataBase (BGRDB, Diehl et al. 2018), as well as Cyg X-1 and Cyg X-3. We use the OSA spimodfit (Strong et al. 2005) to perform these fits.

Given the large unknown in the SPI data reduction processes, and although Cyg X-1 is not formally a source sufficiently bright to present the caveats described by Roques & Jourdain (2019) in both states, we also extract SPI spectra using different methods (see Appendix B). The SPI background variation timescale depends on energy and energy bin size (Siegert et al. 2019); additionally, the absolute fluxes of the sources may change with time, even when they are in the same state. To determine the optimal variation timescales for both background and source at the same time, we therefore perform a grid search. To do so, we define 13 variability timescales, between one pointing (∼0.6 h) and no variability (constant background and source), and double the times in each step (in total 169 grid points). The smaller the variability timescale, the more parameters are included in the fit. Although this results in a better likelihood, it tends to over-sample the data. In order to penalize large numbers of parameters, we search for the minimum Akaike information criterion (AIC Akaike 1974; Cavanaugh et al. 1997) instead, which defines the fit quality of models relative to each other (see the Appendix for examples of AIC maps). From this analysis, we find that in the HSS the background variability timescale is rather constant as a function of energy, changing between 12 h and 72 h (∼1/6–1 INTEGRAL orbit). In the LHS, background variations are more frequent. The background variability timescale is found to decrease down to the pointing timescale below 150–200 keV. However, at these energies the background and source patterns in addition to the large number of fitted parameters are too degenerate to determine a robust estimate for a temporal variability timescale of Cyg X-1.

It appears that switching off a strong source variability provides similar flux normalizations to IBIS/ISGRI spectral extractions (c_ISGRI/c_SPI ≈ 1.0, see Table 2) in the LHS, but to a lesser extent in the HSS (c_ISGRI/c_SPI ≈ 1.3). Above 300 keV (400 keV) in soft (hard) state, the number of fitted parameters for background and source are considerably smaller, also reducing the degeneracy in the maximum likelihood fits. The source spectra are then found to be steady, and the high-energy tail (Sect. 3.2) is constant in time.

We compare our extracted INTEGRAL/SPI hard state spectrum with the INTEGRAL/IBIS Compton-mode spectrum published by R15 in Fig. 3. At the time, the spectrum was blamed for being very hard. Here, we observe a softer spectrum as extracted with SPI compared to the Compton-mode spectrum previously published by R15.

Fig. 3. Comparison of our extracted hard-state INTEGRAL/SPI spectrum with INTEGRAL/IBIS/PICsIT spectra (a) extracted for the same time period as SPI and (b) for data before 2012 as published by Rodriguez et al. (2015).

2.4. Spectral extraction: INTEGRAL/JEM-X

For the Joint European X-ray Monitor spectral extraction, we restrict our selection to scws where the source is in the field of view of the JEM-X unit 1, which is where the source is less than 5° off-axis. This selection results in 836 and 1112 scws for the LHS and HSS, respectively. The data from JEM-X 1 are reduced with OSA v11.0. We follow the standard steps described in the JEM-X user manual⁵. Spectra are extracted in each scw where the source is automatically detected by the software at the image creation and fitting step, and are then computed from ∼3 to ∼34 keV with 32 logarithmic bins. The individual spectra are then combined with the OSAspe_pick tool. The appropriate ancillary response files (arfs) are produced during the spectral extraction and combined with spe_pick, while the redistribution matrix file (rmf) is rebinned from the instrument characteristic standard rmf with j_rebin_rmf. We add 3% systematic error on all spectral channels for each of the stacked spectra, as recommended in the JEM-X user manual, and the spectra are considered between 3 and 25 keV.

3. Results

3.1. Methodology for the spectral analysis

We use both the Interactive Spectral Interpretation System (ISIS; Houck & Denicola 2000) and XSPEC V.12.10 (Arnaud 1996) packages for the spectral analysis.

In order to assess the presence, significance, and main properties of the high-energy tail, we first analyze the data with a simple phenomenological model, that is, a thermal Comptonization continuum in combination with a power law extending into the MeV range (Sect. 3.2). Keeping in mind the suggested origin for the LHS high-energy tail proposed by L11, Jourdain et al. (2012), and R15, we then investigate other possible origins for this component. We use physical models considering the broad-band source spectrum as originating from a single medium composed of a hybrid thermal/nonthermal particle distribution (Sects. 3.3 and 3.4).

3.2. Phenomenological approach: a high-energy tail in both states

The LHS and HSS stacked spectra of JEM-X 1, IBIS/ISGRI, and SPI are shown in Fig. 4. Following the approach presented in J14 and R15, we first focus on the 3–400 keV spectral band. In that energy range, we apply the thermal Comptonization model compTT (Titarchuk 1994) being reflected off the disk, which in xspec notation is

const∗reflect(compTT),

Fig. 4. Long-term (2003–2018), broad-band (3 keV to 2 MeV), and state-resolved (red: HSS, blue: LHS) stacked spectra of JEM-X 1 (light colors), ISGRI (intermediate colors), and SPI (darker colors).

and has been successfully used by J14 and R15. The const component is a multiplicative constant necessitated by the differences in absolute flux calibration between the different instruments used. In the LHS, we freeze the value of the seed photon temperature of the comptt component to kT₀ = 0.2 keV (e.g., Balucinska & Hasinger 1991; Wilms et al. 2006, R15). In the HSS, we also need to add the thermal emission from the disk (diskbb) and an iron line at 6.4 keV (fitted by a gaussian, gauss). Thus, the model is written in xspec as

const∗(reflect(compTT)+diskbb+gauss).

As the spectral resolution of JEM-X prevents a fine description of the iron line, we fix its energy and its width to E_Fe = 6.4keV and σ_Fe = 0.4 keV, respectively. The seed photon temperature of the comptt component is fixed to the temperature of the diskbb component.

Once an appropriate fit has been obtained with the spectra below 400 keV, we then add ISGRI and SPI data above 400 keV and let the parameters vary freely. A clear high-energy excess is visible in both states above the continuum described with the above-mentioned models. We model this additional emission with a power-law component. We assess the significance of the added power-law tail using an F-test. A low P-value probability of the F-test suggests that the null-hypothesis, that is, the assumption that the thermal Comptonization continuum describes the broad-band spectrum well, is false and means that the fit requires the additional power-law component to model the high-energy tail. We find an F-test probability of P = 3.5 × 10⁻²⁴ and 2.6 × 10⁻⁴ in the LHS and HSS, respectively. Table 2 lists all best-fit parameters in the 3–2000 keV range and Fig. 5 shows the best-fit models and their residuals for both states.

Fig. 5. Best-fit models obtained from our phenomenological fitting of the 3–2000 keV LHS (blue) and HSS (red) spectra. JEM-X, ISGRI, and SPI data are represented with light, intermediate, and darker colors respectively. The different model components are shown with different line styles: Comptonized continuum (densely dashed), additional power law (dotted), disk emission (dotted-dashed), and the iron line (loosely dashed). The bottom pan els show the residuals between these models and spectra for the source in both states.

3.2.1. The LHS

In the LHS, the electron temperature, kT ∼ 53 keV, is consistent with the value obtained by J14 and roughly consistent with that of R15. The optical depth, τ ∼ 1.5, seems here slightly higher than in the two previous studies. The reflected photon fraction of Ω/2π ∼ 0.8 is consistent with the one found by J14 but remains higher than in R15 (where Ω/2π = 0.13 ± 0.02). These small differences are probably due to the intrinsically variable nature of the source, our largely extended observations sample and integrated exposure time in this state compared to these previous studies, and the extended spectral range we consider here compared to J14 and R15, and are therefore not surprising. As expected from the comparison of our SPI spectrum with that of R15 (Fig. 3), we find a softer photon index (Γ_tail ∼ 2.06; Table 2) for the additional power-law component than the (very hard) value found by R15 ( = ). Our results nevertheless remain rather consistent with those of J14 ( = 1.8 ± 0.2).

3.2.2. The hard-soft state

Significant residuals arise in the HSS between tens to hundreds of keV for IBIS/ISGRI when extracted with OSA v10.2. We therefore exclude the OSA v10.2 HSS data from the fit. The resulting reduction in count statistics has only a minor impact due to the coverage of the same energy range by SPI.

We are able to integrate a sufficient number of counts in the HSS of Cyg X-1 and can confirm the suggestion previously made by J14 of a possible high-energy tail in the HSS. We find a photon index for the additional power law Γ_tail ≃ 2.0 (Table 2), which is consistent with the value we find in the LHS. We find an electron temperature of , which is not well constrained, and an optical depth of τ ≃ 0.29 consistent with the previous results of J14. We can constrain the reflection component to Ω/2π ≃ 0.5 (Table 2).

3.3. State-resolved spectral modeling: unmagnetized hybrid plasma

3.3.1. Basics of the eqpair model

To go beyond the simple phenomenological approach, we model the data with the emission from a single region, where the electron distribution is hybrid. In such a model, the high-energy tail is explained by the emission of a nonthermal population of accelerated particles. These particles cool down and eventually thermalize because of several processes such as Coulomb collisions. The steady-state distribution is then a hybrid distribution constituted by a low-energy, thermal population producing the bulk of the emission below 400 keV and a high-energy, nonthermal population producing the high-energy tail above 400 keV.

To model such a hybrid plasma, we first use eqpair (Coppi 1999), which has been successfully used in describing the broad-band spectra of Cyg X-1 in different states (Gierliński et al. 1999; McConnell et al. 2002; Wilms et al. 2006; Cadolle Bel et al. 2006; Fritz et al. 2007; Nowak et al. 2011; Del Santo et al. 2013). This model takes into account all microphysical interactions and processes between particles and photons (bremsstrahlung, Compton-scattering, pair production, and electron-electron Coulomb collisions) in a homogeneous, isotropic, and ionized medium of total optical depth τ_T which is the sum of the proton Thomson optical depth τ_p plus the optical depth of electron-positron pairs. Soft photons are injected into the medium as a multicolor disk black-body characterized by the inner disk temperature kT_bb and the total power L_s. Energy is provided to the electron population through two channels: either directly in the thermal pool of electrons with power L_th, or as nonthermal particles with a power-law distribution with index Γ_inj and power L_nth. The total power injected into electrons is L_h = L_th + L_nth. In steady state, the source luminosity is equal to the total injected power: L = L_s + L_h. In the code, the luminosity is described by the compactness parameter:

(1)

where σ_T is the Thomson cross-section, m_e the electron mass, and R is the typical source size. Other compactness parameters l_s, l_h = l_th + l_nth are defined in a similar manner. The overall spectral shape is strongly dependent on the dominance of high-energy emission, that is, the ratios l_h/l_s and l_nth/l_h.

3.3.2. Application to our data

As in all our different approaches, we freeze the system inclination to 27°. In the LHS, we fix the temperature of the disk photons to kT_bb = 0.2 keV and the disk compactness to l_s = 1 (e.g., Gierliński et al. 1999; McConnell et al. 2002; Del Santo et al. 2013). In the HSS, the spectrum is dominated by the emission of the disk photons and therefore, we let the disk temperature, kT_bb, vary freely and fix l_s to 10 (Gierliński et al. 1999; McConnell et al. 2002; Del Santo et al. 2013). We further assume the Lorentz factor of the nonthermal plasma to be power-law distributed between γ_min = 1.3 and γ_max = 1000. The radius of the corona is set to 10⁷ cm and we assume solar abundances ANGR (Anders & Tobin 1989) for all metals in the reflection modeling. The ionization parameter is set to zero and we do not inject pairs. As in our previous approach, we also add a Gaussian to the HSS model.

We show the best-fit model in Fig. 6 and list all parameters in Table 3. Our parameters are consistent overall with the previous studies indicated above.

Fig. 6. Best models obtained from the fits to spectra for Cyg X-1 with the unmagnetized (left) and magnetized (right) hybrid corona models. Plain lines correspond to the full models whereas the loosely dashed line shows the iron line component. JEM-X, ISGRI, and SPI data are represented in light, intermediate, and darker colors, respectively.

The LHS is clearly dominated by Comptonization at an optical depth of ∼1.3 and with a high fraction of accelerated nonthermal particles with l_nth/l_h ∼ 0.9. This results in a spectrum that is dominated by the hard X-rays, that is, l_h/l_s ∼ 6.3. The accelerated electron distribution is a power law with an injection electron index Γ_inj ≃ 2.9 (Table 3).

Our HSS spectra show some emission up to 600 keV, i.e., to slightly higher energies than in most previous studies. In this state, our fit becomes unstable when letting τ_p vary freely. We therefore freeze this parameter to τ_p= 0.2 in line with the values previously obtained with eqpair (Gierliński et al. 1999; Frontera et al. 2001; McConnell et al. 2002; Sabatini et al. 2013; Del Santo et al. 2013) and in coherence with our phenomenological approach (Sect. 3.2.2; Table 2). We note that we need an injection electron index of Γ_inj ∼ 3.0 consistent with the value we find in the LHS (Table 3). The spectrum is dominated by the disk thermal emission as indicated by a l_h/l_s ∼ 0.4 fraction weaker and a photon disk temperature of kT_bb ∼ 0.6 keV, which is higher than in the LHS and is consistent with the results of Gierliński et al. (1999), Frontera et al. (2001), McConnell et al. (2002), Sabatini et al. (2013), Del Santo et al. (2013). The fraction of particles emitted by nonthermal processes l_nth/l_h is consistent with a purely nonthermal distribution (> 77%). Finally, the reflected photon fraction is consistent with our previous approach, although it is less constrained.

3.4. State-resolved spectral modeling: magnetized hybrid plasma

3.4.1. Basics of the belm model

In the previous section, we model the emission from an unmagnetized medium with the eqpair model. It is nevertheless believed that the disk and the jet are significantly magnetized (e.g., Blandford & Znajek 1977; Ferreira et al. 2006). Moreover, nonthermal particles are often assumed to be accelerated by magnetic reconnection, which also requires strong magnetic fields (e.g., Poutanen & Vurm 2009).

In this section we therefore extend the previous modeling by additionally considering magnetic fields as a natural explanation for the accelerated, nonthermal particle distribution that we used for eqpair. To do so, we use the belm code (Belmont et al. 2008), which uses the same parameters as eqpair, but additionally introduces the magnetic field intensity through the magnetic compactness,

(2)

The inclusion of the magnetic field yields several effects over the modeling presented so far. Firstly, additional synchrotron radiation is produced, which can then be Comptonized (synchrotron self-Compton). This directly impacts the spectral shape by adding a new broad-band component and indirectly by increasing the particle cooling rate. Secondly, low-energy particles can synchrotron-absorb high-energy photons, which heats these particles, resulting in efficient thermalization.

To moderate the model degeneracy, we focus on a purely nonthermal model (l_th = 0, see also Del Santo et al. 2013). Thus, the total radiation compactness l = l_s + l_th + l_nth becomes l = l_s + l_nth, and l_h = l_nth.

The table we use is computed with five model parameters: l_nth, l_B/l_nth, τ_p, Γ_inj and kT_in. As shown in Malzac & Belmont (2009), the shape of the spectrum depends on ratios such as l_B/l_nth or l_nth/l and not on the value of the total compactness l. Therefore, its value is fixed to l = 10 in the table which represents the typical value of 0.01L_Edd for both states, while l_s is constrained by enforcing l_s + l_nth = 10. The exact magnetic field compactness depends on the normalization N of the model,

(3)

The magnetic field intensity can then be estimated assuming the distance to the source, D = 1.86 kpc, and the size of the corona, R = 20 R_g. As previously, we fix the value of the injected seed photons to kT_in = 0.2 keV in the LHS and τ_p = 0.2 in the HSS. An iron line fixed at 6.4 keV is also added in the HSS. Here, we describe reflection using the convolution code reflect. In eqpair, the reflection model ireflect is directly implemented as part of the model, with the additional capability to model an ionized medium. In our case, we choose a neutral medium in order to guarantee comparability between belm and eqpair.

3.4.2. Application to our data

Table 4 shows the best-fitting parameters for our modelings with belm. The best-fit models and the residuals between those and the source spectra are shown in Fig. 6.

The LHS is only compatible with weak values of the magnetic compactness parameter (l_B/l_nth < 0.04). This results in a coronal magnetic field B < 4.83 × 10⁵(20R_g/R) G compatible with the results obtained by Malzac & Belmont (2009), Del Santo et al. (2013), and Poutanen & Vurm (2009). Larger fields produce large amounts of soft synchrotron photons that Compton cool the thermal electron population. The parameters obtained in this state (Table 4) are consistent with the parameters found with the unmagnetized model eqpair (Table 3), except for the optical depth for which we observe a larger value with belm. The difference in optical depths may arise from a different treatment of the coronal geometry: belm uses an escape probability corresponding to a pure spherical geometry, while eqpair uses a combination of factors related to spherical and slab geometries.

In the HSS, we also find parameters consistent with the unmagnetized model. The parameters are compatible with the previous results of Malzac & Belmont (2009), Poutanen & Vurm (2009), Del Santo et al. (2013). We note that l_B and l_s both contribute to the corona cooling. l_s is by definition the total amount of cold photons entering in the corona, and thus, the fraction of disk photons that participate in the cooling is 100%. We also test the possibility that a fraction of the disk photons will not interact with the corona by adding a disk black body (diskbb) in the model and investigate whether it has an effect on the value of l_B/l_nth and/or l_nth. The seed photon temperature of diskbb is fixed to the temperature kT_in in the belm model. We find that the uncertainties were slightly different: l_B/l_nth = 1.6 ± 0.4 and with the addition of diskbb compared to and without (Table 4). This latest value on l_B/l_nth leads to a magnetic field of ) G.

The shape of the Comptonization spectrum depends on the heating/cooling ratio; a smaller l_s and a larger l_B can lead to the same Comptonization spectrum. However, contrary to a disk, the thin synchrotron spectrum extends above 1 keV; in the LHS, this component can be hidden by the Comptonized component, but in the HSS, if there was a synchrotron emission comparable to the disk emission, we should see it in the spectrum, which is not the case and this is confirmed by the fits with an additional disk black-body component. We discuss these results and their implications in the following section.

4. Discussion

4.1. Summary of the main results

The extended INTEGRAL data set considered here allow us to confirm the detection of a high-energy tail in the LHS (e.g., Cadolle Bel et al. 2006; L11; J14; R15). We also report the detection of a high-energy tail in the HSS with INTEGRAL, extending to at least 600 keV as seen with SPI (Fig. 4).

Our empirical model (Sect. 3.2) consists of purely thermal Comptonization (compTT) with reflection off neutral gas and an independent power-law component to model any high-energy tail emission. We bring solid constraints for the high-energy tails in both states: with Γ^LHS ≃ 2.03 and Γ^HSS ≃ 2.0 (Table 2) they, in particular, appear to have similar slopes. Their 200–1000 keV fluxes are 3.5 × 10⁻⁹ and 7.1 × 10⁻¹⁰ ergs/cm²/s which are somewhat different.

In an attempt to gain a more physical description of the data, we tested hybrid Comptonization models for an unmagnetized (eqpair) and magnetized (belm) corona. Our results are consistent with those of previous studies (Gierliński et al. 1999; McConnell et al. 2002; Del Santo et al. 2013). The fits with eqpair lead to a higher value of l_h/l_s in the LHS than in the HSS, which indicates that the source emission is dominated by Comptonization in the LHS. Both states are well modeled and we obtain electron injection indices of and . These results are globally consistent with the parameters we obtain with belm. Moreover, belm provides us constraints on the magnetic field. We found an upper limit B < 4.83 × 10⁵(20R_g/R) Gauss in the LHS, and we are able to constrain Gauss in the HSS. Those values are consistent with upper limits found by Malzac & Belmont (2009), Poutanen & Vurm (2009), Del Santo et al. (2013).

4.2. Interpretation: origin of the high-energy emission

The detection of a high-energy tail in both states leads us to consider two different scenarios:

1. The high-energy tails in the LHS and HSS originate from the same medium.

2. The high-energy tails in the LHS and HSS are the signatures of different physical mechanisms, and thus have distinct origins, despite having photon indices of similar values.

4.2.1. Same origin and spectral interpretation

Our results with the semi-phenomenological approach are consistent with what has been previously reported for Cyg X-1: a significantly higher optical depth in the LHS compared to the HSS, which could suggest that the corona shrinks as the system transits from the LHS to HSS (e.g., Zdziarski et al. 1998; Wilms et al. 2006; Del Santo et al. 2013), and an electron temperature around 50 keV in the LHS. However, this approach does not allow us to treat the broadband spectra in a self-consistent way, and prevents us from verifying whether or not a single medium (a corona) can be responsible for the 100–1000 keV emission.

With the physical, broad-band hybrid spectral models eqpair and belm, in contrast, we can study the evolutionary behavior of a corona as a single self-consistent medium. The first conclusion we can draw is that our spectra are well modelled with both the unmagnetized and magnetized models eqpair and belm. In these frameworks, the high-energy emission is produced by a nonthermal electron population in either a nonmagnetized or a magnetized corona in both states.

As belm furnishes precious constraints on the magnetic field in both states, we can, moreover, draw another conclusion which appears intriguing. The ratio between the equipartition magnetic field compactness l_BR and the magnetic field compactness l_B is given by Lightman & Zdziarski (1987):

(4)

We find and for the LHS and HSS, respectively. In the LHS, this low value argues against magnetic reconnection powering the corona, which is consistent with the conclusions of Malzac & Belmont (2009), Poutanen & Vurm (2009), and Del Santo et al. (2013). If electron acceleration is not possible via magnetic reconnection, it thus implies that the high-energy emission must be produced in a second, independent zone, which is not the corona, and could for example be the jets.

In the framework of belm as a one-zone model, the MeV tail is mainly due to Comptonization and would be inconsistent with the observed high degree of polarization claimed for the LHS so far (L11; Jourdain et al. 2012, R15). In the HSS, the magnetic field is large enough to be able to produce magnetic reconnection, and this thus could be the mechanism responsible for the electron acceleration in the corona (Malzac & Belmont 2009; Poutanen & Vurm 2009; Del Santo et al. 2013). To summarize, an unmagnetized corona could be at the origin of the high-energy emission in both states and is thus a plausible scenario for both, while a magnetized corona seems to be working only in the HSS.

4.2.2. Different origin and comparison with previous results

As indicated, the tails in both long-lasting spectral states do not necessarily have the same origin. Evolved radio jets (e.g., Rushton et al. 2012, and references therein), the high polarization fraction of the high-energy tail above 400 keV (L11; Jourdain et al. 2012; R15), and broadband high-energy spectral features associated with the jets of Cyg X-1 (L11; Jourdain et al. 2012; Zdziarski et al. 2012; R15; Walter & Xu 2017) seem to point towards a synchrotron jet emission in the LHS. Several studies investigate this possibility. Zdziarski et al. (2012) apply their jet model and find that such emission requires an electron acceleration at a rather hard power-law index (1.3–1.6), which is difficult to reconcile with standard acceleration mechanisms. Similar conclusions were found in an another source observed with INTEGRAL, namely GRS 1716–246, for which Bassi et al. (2020) used a jet internal shock emission model (Malzac 2013) in order to model the high-energy emission of this source. However, another study on Cyg X-1 shows that the emission during the LHS can be explained by a self-consistent jet model (Kantzas et al. 2020). In their model, a broadband jet-synchrotron component is the origin of both the radio continuum and the MeV tail in the LHS. A photon index of Γ ∼ 2 in the MeV range would be compatible with their results. The jet origin is further supported by the detection between 60 MeV and 20 GeV by Fermi/LAT (Bodaghee et al. 2013; Malyshev et al. 2013; Zanin et al. 2016; Zdziarski et al. 2017), which may also hint towards a hadronic contribution to the jet-accelerated plasma (Kantzas et al. 2020).

A scenario that lies beyond the scope of our paper, and is hard to constrain and fully discuss with our data, is the putative presence of a jet in both states. This possibility arises from recent claims of possible detected jets in the HSS. In particular, shortly after the LHS-to-HSS state transition around MJS 55400, Rushton et al. (2012) detect a weak and compact jet on VLBI scales. It is unclear whether this feature has remained undetected in the long-term stable HSS before, or whether this is a remnant feature of the state transition that happened close-by in time. A recent study claims the detection of radio jets in the HSS of Cyg X-1 (Zdziarski et al. 2020). The authors observe a clear correlation between 15 GHz radio flux and the simultaneous 15–50 keV flux not only in the LHS but also, albeit with a different slope, in the HSS. However, as opposed to the LHS, the soft X-rays at a few keV are reported to be time-lagged by about 100 days in the HSS, suggesting advection from the donor through the accretion disk. The absence of a time-lag between the radio and hard (15–50 keV) X-ray emission in both states might, here, imply a jet origin for the latter hard X-ray radiation. However, we note that the significantly lower degree of polarization below ∼400 keV (Jourdain et al. 2012) and both our spectral modeling and that of Kantzas et al. (2020) suggest that the persistent jet may only contribute above hundreds of keV during the LHS.

No radio emission has been detected in any other microquasar in the HSS down to a very low flux level (e.g., Fender et al. 1999; Remillard & McClintock 2006; Belloni 2010). Still, a recent study suggests that jets could survive in the HSS with heavily suppressed radio emission that may or may not be detectable depending on the observing technique and intrinsic radio-flux level. These so-called “dark jets” are kinematically dominated and radiatively inefficient (Drappeau et al. 2017). In the context of the internal shock model of Drappeau et al. (2017), the radio flux is quenched, but a kinematic jet driven by the disk rather than the corona still subsists. However, this study is based on a broadband modeling of GX 339−4 and H 1743−322, both of which are microquasars with low-mass companion stars and are purely Roche lobe accreting. It is unclear whether or not the results from Drappeau et al. (2017) can be applied to Cyg X-1, a persistently wind-accreting system with a high-mass companion.

5. Conclusion

In this paper, we advance the long-lasting investigation of the nature of the hard MeV tail of Cyg X-1 for both major spectral states. We accumulate data from 15 years of INTEGRAL (JEM-X, IBIS and SPI) observations of the source and reach unprecedented count statistics in the 3–2000 keV range, not only in the LHS, but also in the HSS, which allows us to firmly confirm the suspected presence of a high-energy component —extending to at least 600 keV– in the spectra of the source.

We follow two approaches with the modeling of these data. In the first, empirical model, we fit a thermal Comptonization spectrum independent of a simple power law, describing the MeV tail between hundreds of keV and ∼1 MeV. The slopes of the MeV tails appear compatible between the LHS and the HSS. The origin of this high-energy emission can be tackled with single physical spectral models where the data are self-consistently considered. We therefore used the hybrid thermal/nonthermal Comptonization codes eqpair and belm, accounting for an unmagnetized and magnetized single-zone corona, respectively. We reasonably successfully reproduce the spectral changes across the states of Cyg X-1 overall. Within this single-zone approach, an unmagnetized corona is able to account for the high-energy emission in both states. Considering the magnetized coronal model, in the LHS, the magnetic field is too weak to be able to accelerate electrons via magnetic reconnection. In the HSS, this mechanism could be responsible for the nonthermal emission we observe (Malzac & Belmont 2009; Poutanen & Vurm 2009; Del Santo et al. 2013).

We then also discuss our results in the context of physical broadband multi-zone models from the literature and their ability to reproduce the radio-to-gamma-ray emission of Cyg X-1. In the LHS, a persistent radio jet can explain both the radio emission and the MeV tail via a direct synchrotron component, while the X-rays at intermediate energies are well described with hybrid Comptonization. When switching to the HSS, we suggest that the persistent jet may transition to a quenched or failed jet, leaving behind a hybrid corona that could either be magnetized or not, and which is responsible of the nonthermal emission.

Our study does not provide conclusive evidence in favor of any particular one of the scenarios proposed to explain the origin of the high-energy tail in the HSS. A polarimetric analysis of the MeV tails for both states is in preparation, extending the work previously presented in R15. Polarimetry will be key to breaking degeneracies arising between the different physical models, but is currently still affected by large uncertainties. Continued observations and further studies of Cyg X-1, especially in the HSS at energies beyond 1 MeV with INTEGRAL, will be crucial to increasing the count statistics and thus to better constraining the slope of the MeV tail. Complementary studies of other sources would facilitate our understanding of universal properties and differences in such high-energy tails beyond Cyg X-1.

Acknowledgments

We thank D. Kantzas, M. Lucchini, A. Chhotray, and S. Markoff for extensive discussions on the origin of the high-energy tail emission based on physical jet models developed at the Anton Pannekoek Institute for Astronomy, Univ. of Amsterdam. F. C., J. R., P. L. acknowledge partial funding from the French Space Agency (CNES), and the French programme national des hautes énergies (PNHE). F. C. and V. G. acknowledge support from the ESTEC Faculty Visiting Scientist Programme. TS is supported by the German Research Society (DFG-Forschungsstipendium SI 2502/1-1). VG is supported through the Margarete von Wrangell fellowship by the ESF and the Ministry of Science, Research and the Arts Baden-Württemberg. This work was partly funded by Deutsches Zentrum für Luft- und Raumfahrt grant 50 OR 1411 and 50 OR 1909. We made use of ISIS functions provided by ECAP/Remeis Observatory and MIT (http://www.sternwarte.uni-erlangen.de/isis/). We thank J.E. Davis for the development of the slxfig module that has been used to prepare the figures in this work. This work is based on observations with INTEGRAL, an ESA project with instruments and science data centre funded by ESA member states (especially the PI countries: Denmark, France, Germany, Italy, Switzerland, Spain), and with the participation of the Russian Federation and the USA. This research has made use of NASA’s Astrophysics Data System Bibliographic Services (ADS).

References

Appendix A: Variability timescales in Cyg X-1: the AIC map (SPI)

Figure A.1 shows AIC values as a result of maximum likelihood fits with different numbers of parameters for the background variability (timescale: the more parameters, the more variable) and the source (Cyg X-1) variability timescale, which was carried out in order to determine the flux and its uncertainty in each bin. In the HSS, the AIC maps appear with a clear minimum (most probable solution according to AIC), slowly changing the source variability from ∼10 h to constant as a function of energy. In the LHS, the AIC minimum is found at the smallest timescale possible (one scw) for both source and background. For higher energies, the background still appears to vary on the shortest possible timescale, and the source seems constant. This appears unreasonable as the background variability should not change as a function of energy. We believe that the LHS spectral extraction suffers from a degeneracy in the determination of the number of required parameters. We therefore “calibrate” our spectrum using the IBIS spectral extraction for the LHS: using a source model without temporal variability but allowing the background to freely vary leads to an absolute normalization and shape close to the IBIS extraction.

Fig. A.1. Example of AIC maps with a different number of parameters for the background and the source variabilities.

Appendix B: Comparison of different SPI spectral extraction methods

Here, we present two other spectral extraction methods for the SPI data.

B.1. Method 2: Using different event selection

In this method, we consider the threshold of the PSD events, as set in the SPI telescope electronics. These are from 400 keV to 2.2 MeV after revolution 1170 and from nearly 650 keV to 2.2 MeV before. In this energy range, we used the PSD events only. This selection has been prescribed in Roques & Jourdain (2019) in order to avoid spectral artifacts caused by potential spurious events. These spurious events appear at high energies (> 400 keV) for very bright sources (typically above 1 Crab). They are caused by a displacement of low-energy photons toward high energy, when the source is so bright that many photons can be detected in the same detection electronics time window. The energy reconstruction for the SE chain is then not reliable. This does not happen in the PE chain, which records the same event much more rapidly.

B.2. Method 3: Using SPIDAI

In this method, we use the SPI Data Analysis Interface⁶ (SPIDAI) to reduce the data from SPI. We created a sky model containing Cyg X-1 and Cygnus X-3 and set the source variability to ten scws for Cygnus X-3 and five scws for Cyg X-1 (e.g., Bouchet et al. 2003; Cangemi et al. 2021). We then created the background model by setting the variability timescale of the normalization of the background pattern to five scws. In order to avoid the effects due to solar flares, radiation belt entries, and other nonthermal incidents, we removed scws for which the reconstructed counts compared to the detector counts give a χ² > 1.5. This selection reduces the total number of scws by ∼10%. The shadowgrams were then deconvolved to obtain the source flux, and spectra were extracted between 20 and 2000 keV using 27 logaritmically spaced channels.

Given the large number of scws in both states, we extracted six and four different spectra for the LHS and HSS, respectively, each of them was extracted from data of 500–600 scws and corresponds to a distinct period. For each spectrum, we applied a correcting factor of 1/0.85 [above 400 keV or 650 keV according to the considered period, Roques & Jourdain 2019) to account for the change efficiency. We then computed two average spectra for both states. In the HSS, the normalization of the spectrum has a slight variability over the four different periods which explain the larger error bars that we observe in this state when computing the average spectrum.

B.3. Comparison between the methods

The right panel on Fig. B.1 shows the broad SPI spectra (20–2000 keV) used in this paper along with the spectra extracted with the aforementioned methods. The left panel on Fig. B.1 shows the same spectra between 400 keV and 700 keV, that is to say the range for which the selection between SEs and PSD events is different from the the method described in this paper and method 2, described in Appendix B.1.

Fig. B.1. Left panel: spectra extracted with three different methods: method 1 described in Sect. 2.3 of this paper (dark red), method 2 described in Appendix B.1 (medium red), and method 3 using SPIDAI (light red). Right panel: focus on the 400–700 keV energy range.

In Fig. B.1, we show the residuals we obtain by fitting the spectra with a simple cutoff powerlaw (constant*cutoffpl, , ). By fixing the constant parameter of the spectra used in this paper to 1, we obtain constant values of c_method 2 = 0.97 ± 0.02 and c_method 3 = 1.06 ± 0.03 for the LHS and c_method 2 = 1.00 ± 0.02 and c_method 3 = 0.86 ± 0.06 for the HSS. For both states, the three spectra are consistent and we do not observe any significant differences, in particular in the 400–650 keV range, the range most affected by spectral artifacts caused by potential spurious events.

We also investigate whether the parameters extracted in this study would be different using method 2 or method 3 spectrum. We find that the parameters remain unchanged and zre consistent within the uncertainties for both states when applying the semi-physical model or the hybrid corona models (eqpair and belm).

All Tables

All Figures

Fig. 1. Spectral state evolution of Cyg X-1 between 2003 and 2017 after Grinberg et al. (2013). The HSS is shown in red, the IS in green, and the LHS in blue. Unclassified data are plotted in gray. Upper panel: 30–50 keV INTEGRAL/ISGRI light curve. The vertical dashed line indicates the change of the INTEGRAL data analysis software. Lower three panels: RXTE/ASM, MAXI, and Swift/BAT light curves of Cyg X-1.

In the text

	Fig. 2. SPI exposure of Cyg X-1 in its soft state in units of 106 cm2 s−1. The fully coded field of view of SPI is 16° ×16°. Green points indicate point sources that have to be taken into account in the SPI data analysis. Red points show other sources listed in the SPI catalog of Bouchet et al. (2011) but only seen below ≲100 keV or at specific flare times.
In the text

	Fig. 3. Comparison of our extracted hard-state INTEGRAL/SPI spectrum with INTEGRAL/IBIS/PICsIT spectra (a) extracted for the same time period as SPI and (b) for data before 2012 as published by Rodriguez et al. (2015).
In the text

	Fig. 4. Long-term (2003–2018), broad-band (3 keV to 2 MeV), and state-resolved (red: HSS, blue: LHS) stacked spectra of JEM-X 1 (light colors), ISGRI (intermediate colors), and SPI (darker colors).
In the text

Fig. 5. Best-fit models obtained from our phenomenological fitting of the 3–2000 keV LHS (blue) and HSS (red) spectra. JEM-X, ISGRI, and SPI data are represented with light, intermediate, and darker colors respectively. The different model components are shown with different line styles: Comptonized continuum (densely dashed), additional power law (dotted), disk emission (dotted-dashed), and the iron line (loosely dashed). The bottom pan els show the residuals between these models and spectra for the source in both states.

In the text

	Fig. 6. Best models obtained from the fits to spectra for Cyg X-1 with the unmagnetized (left) and magnetized (right) hybrid corona models. Plain lines correspond to the full models whereas the loosely dashed line shows the iron line component. JEM-X, ISGRI, and SPI data are represented in light, intermediate, and darker colors, respectively.
In the text

	Fig. A.1. Example of AIC maps with a different number of parameters for the background and the source variabilities.
In the text

	Fig. B.1. Left panel: spectra extracted with three different methods: method 1 described in Sect. 2.3 of this paper (dark red), method 2 described in Appendix B.1 (medium red), and method 3 using SPIDAI (light red). Right panel: focus on the 400–700 keV energy range.
In the text