Potential origin of the state-dependent high-energy tail in the black hole microquasar Cygnus X-1 as seen with INTEGRAL

0.1-10 MeV observations of the black hole microquasar Cygnus X-1 have shown the presence of a spectral feature in the form of a power law in addition to the standard black body and Comptonization components observed by INTEGRAL. This so-called"high-energy tail"has recently been shown to be strong in its hard spectral state and interpreted as high-energy part of the emission from a compact jet. This result was, however, obtained from a data set dominated by hard state observations. In the soft state, only upper limits on the presence and hence the potential parameters of a high-energy tail could be derived. Using an extended data set we aim at obtaining better constraints on the properties of this spectral component in both states. We make use of data obtained from 15 years of observations with the INTEGRAL satellite. The data set is separated into the different states and we analyse stacked state-resolved spectra obtained from the X-ray monitors, the gamma-ray imager, and the gamma-ray spectrometer onboard. A high-energy component is detected in both states confirming its earlier detection in the hard state and its suspected presence in the soft state with INTEGRAL. We first characterize the high-energy tail components in the two states through a model-independent, phenomenological analysis. We then apply physical models based on hybrid Comptonization. The spectra are well modeled in all cases, with a similar goodness of the fits. While in the phenomenological approach the high-enery tail has similar indices in both states, the fits with the physical models seem to indicate different properties. We discuss the potential origins of the high-energy components in both states, and favor an interpretation where the part of the high-energy component is due to a compact jet in the hard state and hybrid Comptonization in either a magnetised or non-magnetised corona in the soft state.


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 domi-Send offprint requests to: fcangemi@lpnhe.in2p3.fr ; tbeuchert@eso.org nated 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). A&A proofs: manuscript no. Corrected6 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 stateresolved 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 INTE-GRAL (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, Miller-Jones et al., submitted), 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 (Mc-Connell 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 , 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.

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, Fig.  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.  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 −1 for OSA v10.2 and ∼134 counts s −1 for OSA v11. IBIS/ISGRI data were reduced with the standard procedure 3 . 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.

Spectral extraction: INTEGRAL/IBIS/ISGRI
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 3 . 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.

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 timedependent 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. Cy g X− 2 IG R J2 13 35 +5 10 5 Cy g X− 1 Cy g X− 3 We consider single events (SEs) below 530 keV and pulseshape 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 4 . 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 4 .
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.

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 5 . 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 OSA spe_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.

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 (Sect. 3.3 and 3.4).

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), 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 = 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.4 keV 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 highenergy tail. We find an F-test probability of P = 3.5 × 10 −24 and 2.6 × 10 −4 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.

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 (Γ R15 tail =1.4 +0.2 −0.3 ). Our results nevertheless remain rather consistent with those of J14 (Γ J14 tail =1.8 ± 0.2). 10 2 10 3 Energy (keV)

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 62 +71 −28 , 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) .

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: 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 highenergy emission, that is, the ratios l h /l s and l nth /l h .

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;Mc-Connell et al. 2002;Del Santo et al. 2013). We further assume the Lorentz factor of the nonthermal plasma to be power-law dis- tributed between γ min = 1.3 and γ max = 1000. The radius of the corona is set to 10 7 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.
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 (Sec. 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.

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, 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 .
Article number, page 7 of 14 A&A proofs: manuscript no. Corrected6 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, 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. 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 5 (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.

Application to our data
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 l nth = 2.2 +0.9 −0.2 with the addition of diskbb compared to l B /l nth = 1.6 +0.4 −0.7 and l nth = 2.2 +0.7 −0.2 without (Table 4). This latest value on l B /l nth leads to a magnetic field of B = 1.9 +1.6 −1.3 × 10 6 (20R g /R) 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.

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 −9 and 7.1 × 10 −10 ergs/cm 2 /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 Γ LHS inj 2.9 and Γ HSS inj 3.1. 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 5 (20R g /R) Gauss in the LHS, and we are able to constrain B = 1.9 +1.6 −1.3 × 10 6 (20R g /R) Gauss in the HSS. Those values are consistent with upper limits found by Malzac & Belmont (2009), Poutanen & Vurm (2009), Del Santo et al. (2013.

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

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 B R and the magnetic field compactness l B is given by Lightman & Zdziarski (1987): We find l B l B R < 0.08 and 1.4 +0.7 −0.8 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.

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 IN-TEGRAL, 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.

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. Left pannel: Spectra extracted with three different methods: method 1 described in Sec. 2.3 of this paper (dark red), method 2 described in Sec. B.1 (medium red), and method 3 using SPIDAI (light red). Right pannel: Focus on the 400-700 keV energy range.