© ESO 2020

1. Introduction

Cyg X–3, V4641 Sgr, V404 Cyg and GRS 1915+105 are unique sources even in the fairly non-homogenous group of X-ray binaries (XRBs). They present very different companion stars: Cyg X–3 harbors a Wolf-Rayet companion (van Kerkwijk et al. 1996; Koljonen & Maccarone 2017), making it a high-mass XRB, while V4641 Sgr, V404 Cyg, and GRS 1915+105 are low-mass XRBs with a late B-type (Orosz et al. 2001) star, a K-type subgiant (Casares et al. 1992; King 1993; Khargharia et al. 2010), and a K-type giant star (Greiner et al. 2001) as donors, respectively. However, they share some similarities that are unique among the XRB population. They are all very powerful X-ray emitters, with Cyg X–3 persistently emitting a luminosity of 10³⁸ erg s⁻¹ in the X-ray band, V4641 Sgr and V404 Cyg exhibiting luminous outbursts where the X-ray luminosity can exceed the Eddington luminosity for a 10 solar mass black hole (Revnivtsev et al. 2002; Motta et al. 2017a), and GRS 1915+105 being in outburst for the past 27 years (Castro-Tirado et al. 1994; Fender & Belloni 2004), with luminosities reaching and surpassing the Eddington limit (e.g., Done et al. 2004).

Except for Cyg X–3, all these sources have long orbital periods and thus large accretion disks, and there is evidence of a high orbital inclination. Therefore, any large-scale geometrical change in the accretion disk such as puffing-up or warping of the accretion flow, or an equatorial outflow can cause local obscuration events. Evidence of this could be seen in the June 2015 outburst of V404 Cyg, which showed highly variable high column density material absorbing the X-ray continuum that remained hard throughout the outburst (Motta et al. 2017a). Especially so-called X-ray plateaus with diminished X-ray luminosity and softer spectra suggested a heavy obscuration of the intrinsic emission (Motta et al. 2017b; Sánchez-Fernández et al. 2017). Similarly, the high-luminosity active accretion phases of V4641 Sgr can be very rapid with a heavily absorbed hard X-ray continuum (Muñoz-Darias et al. 2018; Revnivtsev et al. 2002). In both systems, there is evidence that the intense likely super-Eddington X-ray emission drives a strong disk wind, thus expelling a significant amount of mass to surround the systems (King et al. 2015; Muñoz-Darias et al. 2018).

On the other hand, Cyg X–3 orbits its companion star at a close distance with a short 4.8 h orbital period (Parsignault et al. 1972). Because the companion is a Wolf-Rayet star exhibiting a heavy stellar wind that extends much farther than the binary orbit, Cyg X–3 is constantly embedded in a high-density environment that affects its X-ray spectra in all accretion states (Szostek & Zdziarski 2008; Zdziarski et al. 2010; Koljonen et al. 2018). This material is optically thick in X-rays as a result of the absorption of metals and Compton scattering, causing iron absorption edges (Koljonen et al. 2018) and Compton downscattering (Zdziarski et al. 2010) and/or Compton scattering out of the line of sight to the intrinsic X-ray continuum. Highly ionized iron lines are resolved with Chandra and reveal a distinct component of gas at much higher ionization, in addition to a component from fluorescence by neutral or near-neutral material. Attempts to unify the iron emission with that of lower-Z elements implied a need for an additional absorption component, possibly associated with a disk wind (Kallman et al. 2019).

GRS 1915+105 is known to have ionized accretion disk wind (Neilsen & Lee 2009), but a high obscuration like that in the other three sources has not been observed. However, GRS 1915+105 recently entered a new accretion state that presents lower fluxes throughout its spectral energy distribution than ever before during its 27 year-long outburst. In this state, sporadic X-ray flares have been observed (e.g., Iwakiri et al. 2019; Neilsen et al. 2019; Jithesh et al. 2019) in addition to the X-ray spectra, indicating heavy obscuration (Miller et al. 2019). Strong radio flares were also observed in the flaring period (Motta et al. 2019; Trushkin et al. 2019; Koljonen et al. 2019), indicating episodic jet emission that is also similar to the multiwavelength evolution of V404 Cyg and Cyg X–3. The similarity of the X-ray spectra in Cyg X–3 and GRS 1915+105 has previously been noted by Vrtilek & Boroson (2013) and Zdziarski et al. (2016), who studied the color-color-intensity diagrams of XRBs and found that GRS 1915+105 and Cyg X–3 occupy an area that is different from that of other black hole or neutron star XRBs. This underlines the connection between the two and their likely “messy” surroundings.

In this paper, we study spectra obtained with the Nuclear Spectroscopic Telescope Array (NuSTAR) and Rossi X-ray Timing Explorer (RXTE) of the four luminous hard X-ray sources Cyg X–3, V404 Cyg, V4641 Sgr, and GRS 1915+105, which are all likely surrounded or occluded in the line of sight by dense material that affects the intrinsic X-ray emission at specific times. The data processing of all sources is described in Sect. 2. In Sect. 3 we show that the X-ray spectra are very similar and peculiar at specific times in the evolution in all four sources, which is caused by the high-density environments in which they are embedded. The mutual spectral characteristics include a low-energy cutoff of the hard X-ray spectra, absorption edges of highly ionized iron, ionized or Doppler-shifted iron emission or absorption lines, and high absorption. We furthermore fit the time-averaged spectra, and also time-resolved spectra in case of V404 Cyg, with two physically motivated models that either describe a scenario in which all the intrinsic emission is reprocessed in the surrounding matter or in which the emitter is surrounded by a thick torus with variable opening angle. This underlines the assumption of X-ray obscuration in all sources. Using the results from both fits, we discuss in Sect. 4 that the (outflowing) obscuring matter in V404 Cyg and GRS 1915+105 shows a change in geometry that is linked to the radio (jet) evolution observed from the sources, in addition to a change in the intrinsic X-ray emission. The sources display X-ray obscuration with varying intrinsic luminosities from lower than 1% of the Eddington luminosity up to the Eddington limit within the framework of the models. This indicates that different factors cause the obscuration. We therefore further discuss the effect of these results for other sources. Finally, we conclude in Sect. 5.

2. Observations and data reduction

2.1. NuSTAR

In the case of V404 Cyg, we selected the NuSTAR observation that was taken during the June outburst of 2015 nearing the end of the 12 day flaring activity that contained both a plateau spectral state with slow spectral evolution and low flux density and a flaring state with rapid spectral changes and highly variable flux density (pointing 90102007002). This observation was previously analyzed in Walton et al. (2017), although they mostly concentrated on analyzing the low-absorption flaring periods with high count rates. For GRS 1915+105, we selected observations that were taken when the source descended to a very anomalous low-flux state in June 2019 that was interspersed by very luminous flares (pointings 90501321002, 30502008002, and 30502008004). For Cyg X–3, we selected the only NuSTAR observation that was taken in the hard state (pointing 10102002002). All observations used in this paper are tabulated in Table 1.

We reduced the NuSTAR data from the two focal plane modules (FPMA and FPMB) using NUPIPELINE. We used a circular source region with a 100 arcsec radius centered on the location of the source, and circular background regions with a 100 arcsec radius that were selected from a sourceless region in the detector image. The source region size was a compromise of including most of the point-spread function but avoiding to confuse it with the possible contribution from the scattering halo that is present in the data of V404 Cyg (Beardmore et al. 2016; Heinz et al. 2016; Vasilopoulos & Petropoulou 2016). The pipeline was run with the parameters TENTACLE = “YES” and SAAMODE = “OPTIMIZED”. The former requires a simultaneous increase in the CdZnTe detector event count rates and the observed shield single rates, and the latter allows identification and flagging of time intervals in which the CdZnTe detector event count rates show an increase when the spacecraft enters the South Atlantic Anomaly (SAA). The data reduction was performed with HEASOFT 6.26.1.

We extracted an averaged spectrum from the two detectors using whole pointings in the case of Cyg X–3 and GRS 1915+105, and spectra with several intervals in the case of V404 Cyg (see below). The broadband (3–79 keV) NuSTAR count rate ranged between 51–188 cts s⁻¹ for Cyg X–3, which mostly arises from the orbital modulation (a factor of 2–3). The total exposure (∼21 ks) is longer than the orbit (∼17.3 ks), therefore we can expect the effect of orbital modulation on the spectral components to average out. For GRS 1915+105, the NuSTAR count rate was found to vary between 1 and 81 cts s⁻¹ during the ∼30–60 ks long pointings, which is considerably lower than usually observed. GRS 1915+105 is famous for its plethora of X-ray variability states (Belloni & Altamirano 2013); the flux and hardness vary on short timescales. For the observations considered here, the hardness ratio between 3–5 keV and 10–79 keV remained relatively constant with 0.8 ± 0.1, 0.22 ± 0.03, and 0.2 ± 0.1 for pointings 90501321002, 30502008002, and 30502008004, respectively. We therefore consider the averaged spectrum as a relatively accurate representation of the spectral shape, although we cannot rule out fast changes in the source spectrum.

For V404 Cyg, we concentrated on analyzing the times preceding and in between the intense X-ray flaring with the total count rate not exceeding ∼1000 cts s⁻¹ in order to study spectral characteristics when the intrinsic X-ray emission was likely obscured (the low-absorption intense flaring spectra were studied in detail in Walton et al. 2017). Because the spectrum changed rapidly during the pointing, we analyzed the spectral evolution in several ways. We divided the pointing into 70 segments, where the extracted spectrum has 30 000 counts to ensure sufficient spectral quality. We excluded spectra from the analysis that contained count rates exceeding 1050 cts s⁻¹ in order to fully include the preflaring period (the first ∼21 ks of the pointing) and periods with low count rates between the luminous flares. This resulted in 25 individual time bins ranging from 420 s to 3980 s in exposure time. We furthermore binned time-resolved spectra that were similar in shape to spectral epochs in order to increase the spectral quality for detailed modeling.

For X-ray modeling, we binned the data to a minimum signal-to-noise ratio (S/N) of 30 in the full band 3–79 keV. The spectral fitting was performed using the Interactive Spectral Interpretation System (ISIS; Houck 2002). In the modeling, a constant factor was added to account for the flux difference of the NuSTAR detectors. In some pointings (e.g., 30502008004 of GRS 1915+105), the discrepancies in the FPMA and FPMB data cannot be explained by a simple constant, especially in the 6–10 keV region. This can affect the fit quality significantly. Therefore we also fit the same models to the spectra from a single detector and present its fit quality as well. However, all the model parameters we present here are estimated from fits to data from both detectors. All the NuSTAR fluxes are normalized to the FPMA detector.

For the X-ray timing analysis, we extracted 2⁻⁶ s light curves from three energy bands: 3–10 keV, 10–79 keV, and 3–79 keV. The cospectra, which can be used as a proxy for white-noise-subtracted power spectral densities (PSDs), were calculated from 512 s long segments averaged over each good time interval (GTI) using Matteo’s Libraries and Tools in Python for NuSTAR timing (MALTPYNT; Bachetti 2015). The cospectrum is used to mitigate instrumental effects in the NuSTAR light curves (Bachetti et al. 2015). We used the rms normalization and binned the cospectra geometrically by a factor of 1.1–1.5 before importing them to ISIS for model fitting.

2.2. RXTE

We downloaded all the proportional counter array (RXTE/PCA) data from the High Energy Astrophysics Science Archive Research Center (HEASARC) during outbursts of V4641 Sgr in 1999, 2002, 2003, and 2005, and selected two representative observations for spectral modeling that resemble the data from Cyg X–3 and V404 Cyg (pointings 80054-08-01-01 and 70119-01-01-14). The former was taken during the outburst of 2003 and was analyzed in Maitra & Bailyn (2006). However, only the first ∼2 ks were studied in their paper, and here we concentrate on the latter part of the light curve. The other selected pointing was taken during the outburst of 2002, and we are not aware that this spectrum has been studied in detail elsewhere. The pointing 80054-08-01-01 was taken at MJD 52857.37 with an exposure of 2.9 ks and a mean count rate of the proportional counter unit 2 (PCU2) of 67 cts s⁻¹, while the pointing 70119-01-01-14 was taken at MJD 52417.81 with an exposure of 0.7 ks and a mean PCU2 countrate of 73 cts s⁻¹. Because neither light curve presented significant spectral changes during the pointings, we extracted the average spectrum for spectral modeling in both cases.

RXTE/PCA data were reduced using the methods described in the RXTE cookbook with HEASOFT 6.26.1. The 128-channel energy spectra were extracted from the standard-2 data using all available PCUs and all layers. For the spectral fitting, we ignore bins below 3 keV and above 40 keV, binned the data to S/N = 5.5, and added 0.5% systematics to each channel. For the timing analysis, we extracted 0.125 s light curves from standard-1 data. We calculated the averaged PSD using 512 s light-curve segments and binned them geometrically by a factor of 1.1 before importing them to ISIS for model fitting.

3. Results

3.1. X-ray spectra: Overview

Figure 1 shows some of the X-ray spectra from the observations tabulated in Table 1. In the top panel, the hard-state spectra of V404 Cyg (from the beginning of the pointing before spectral softening, see Sect. 3.3.2), Cyg X–3, and a 2002 outburst peak spectrum of V4641 Sgr show a striking similarity in spectral shape with a similar absorption profile, a broad peak at the iron line region, a sharp drop at the iron edge energies, and a low-energy cutoff in the hard X-rays. This type of spectrum is not observed from any other XRB. The closest comparison can be found in Compton-thick active galactic nuclei (AGN; e.g., Baloković et al. 2014; Bauer et al. 2015). The GRS 1915+105 spectrum from the anomalous low-luminosity X-ray state before the X-ray or radio flaring shows similar features, but with a prominent iron absorption line and a higher X-ray cutoff energy. The spectra in Fig. 1 (top panel) display different flux densities: the V404 Cyg is a factor of ∼2, ∼6, and ∼20 brighter than Cyg X–3, GRS 1915+105, and V4641 Sgr, respectively. However, when the distances (see Table 2) are taken into account, the X-ray luminosity of V404 Cyg is a factor of ∼3 brighter than V4641 Sgr, but a factor of ∼2 and ∼4 dimmer than GRS 1915+105 and Cyg X–3, respectively, in this state.

Fig. 1. Collection of X-ray spectra from the four sources ordered into two panels according to their spectral shape. Top: NuSTAR spectra (FPMA/FPMB data both included) of V404 Cyg during the outburst period of 2015 (preceding intense X-ray flaring period), Cyg X–3 in the hard state, GRS 1915+105 during the anomalous low-luminosity state of 2019–2020 (preceding a period of sporadic X-ray and radio flaring), and the RXTE/PCA spectrum of V4641 Sgr during the outburst period of 2002. Bottom: NuSTAR spectra of V404 Cyg during the outburst period of 2015 (in between intense X-ray flaring), RXTE/PCA spectrum of V4641 Sgr during the outburst period of 2003, and NuSTAR spectra of GRS 1915+105 from the anomalous low-luminosity state (after the intense X-ray and radio flaring period). Note the spectral similarity between sources. The spectra of GRS 1915+105 and V4641 Sgr are renormalized by the amount shown for illustrative purposes.

In the bottom panel of Fig. 1, a harder X-ray spectrum is shown for V404 Cyg (during X-ray flaring) and V4641 Sgr (2003 outburst peak spectrum), which share approximately the same shape, in addition to a spectrum of GRS 1915+105 after a high-intensity X-ray and radio-flaring period in the low-luminosity state. However, no similar spectrum can be found for Cyg X–3 because all the other spectral states are softer (see all the different spectra from different accretion states in Koljonen et al. 2010). This might therefore indicate that these accretion states are less absorbed or obscured, and the spectral shape might be explained by strong reflection from the accretion disk surface. However, we show here that the observed spectral evolution is also compatible with a change in the geometry of the obscuring material, or with a change in the ionization structure.

All the spectra show a strong iron line with a line width above 200 eV. This has been resolved as a combination of a neutral iron Kα line (6.4 keV) and ionized iron lines, with the strongest component arising from the Fe XXV Kα line (6.7 keV) and possibly from Fe XXVI Lyα line (7.0 keV) with Chandra in the case of Cyg X–3 (Paerels et al. 2000; Kallman et al. 2019) and V404 Cyg (King et al. 2015). The energy band from ∼7 keV to ∼10 keV appears to be affected by the absorption of the above-mentioned species of iron in all sources. Their ionization energies are 7.1 keV, 8.8 keV, and 9.2 keV for Fe Kα, Fe XXV Kα, and Fe XXVI Lyα, respectively. In addition, there is a strong visible iron absorption line around 6.5 keV in the preflare spectrum of GRS 1915+105, and a weaker line in the preflare spectrum of V404 Cyg. Moreover, there might be an indication of emission from an iron Fe XXV Kβ line (7.8 keV) or Ni XXVII Kα (7.8 keV) in the Cyg X–3 spectrum and the Fe XXVI Lyβ line (8.3 keV) in the GRS 1915+105 preflare spectrum. Curiously, some of the V404 Cyg spectra show iron line centroids close to 6.3 keV, similar to findings of Motta et al. (2017a), indicating that the redshifted neutral iron Kα line arises either very close to the compact object (gravitational redshift) or in a medium moving away (Doppler redshift). The strong photoionized emission and absorption lines point to a significant amount of absorbing matter in the line of sight to all sources.

All spectra exhibit a strong curvature or a cutoff in the hard X-rays around 20–30 keV. This is atypical for a hard-state XRB. This cutoff might be caused by a heavy absorption of the intrinsic (cutoff) power-law spectrum, or by a Compton downscattering in the accretion disk (i.e., reflection spectrum) or the surrounding medium. In addition, the spectra shown in Fig. 1 (top panel) present a rather sharp upturn at ∼10 keV, possibly indicating a location where two model components meet that can have very different spectral slopes. In the following, we consider these scenarios by fitting the data with appropriate spectral models.

3.2. X-ray spectra: Initial modeling

In the hard X-ray spectral state, the hard X-ray spectrum of XRBs can typically be fit with a cutoff power law or a Comptonized continuum with a power-law photon index in the range of Γ = 1.5–2.0 and cutoff energy of ∼100 keV. The spectra presented in Fig. 1, especially in the top panel, are quite unlike typical XRB hard-state spectra. Thus, we started by finding and fitting curved models to the hard X-ray data (above 10 keV) of V404 Cyg (preflare or plateau spectrum, as shown in Fig. 1, top panel). The resulting model parameters are tabulated in Table 3 and the models are shown in Fig. 2 (left panel).

Fig. 2. Initial modeling of the V404 Cyg preflare spectrum. Left: best-fit models fit to the hard X-ray data (10–79 keV), but plotted in the full data range. Different models are labeled, and the parameters are tabulated in Table 3. Solid lines refer to pure continuum models, dashed lines to absorbed continuum models, and dotted lines to reprocessed continuum models. Middle: fitting an absorbed cutoff power-law model to the full data range. The dot-dashed blue line corresponds to model A2 fit to the hard X-ray data (see left panel), the dotted green line corresponds to absorption and disk blackbody components added to the model, the dashed yellow line corresponds to partial absorption and smeared edge components added to the model, and the solid red line shows emission and absorption lines added to the model (parameters of the final model are tabulated in Table 4). See the text for more details. Right: fitting a reprocessed thermal Comptonization model to the full data range. The dot-dashed blue line corresponds to model R3 fit to the hard X-ray data (see the left panel), the dotted green line shows an increased value of the ionization parameter, the dashed yellow line corresponds to absorption and smeared edge components added to the model, and the solid red line shows partial absorption and an absorption line added to the model (parameters of the final model are tabulated in Table 4). See the text for more details.

To model the strong spectral curvature, we first tried a cutoff power-law model (C1 in Fig. 2 and Table 3). The best-fit model is highly inverted, with a low-energy cutoff (power-law index, Γ ∼ −1.6, cutoff energy, E_cut ∼ 6 keV), and the fit is poor (${\chi }_{\mathrm{red}}^2=2.6$). Obviously, the power-law index is too low for any reasonable physical scenario. Next, we tried a thermal Comptonization model in a spherical plasma cloud (COMPTT; Titarchuk 1994). The resulting model (C2 in Fig. 2 and Table 3) has an electron temperature of kT_e ∼ 24 keV, a seed photon temperature of kT_s ∼ 5 keV, and an optical depth of τ ∼ 0.7. Model C2 provides a better fit to the data (${\chi }_{\mathrm{red}}^2=1.8$) than model C1, but it is still poor. A tendency for a high seed-photon temperature in fitting thermal Comptonization models to the hard X-ray data of V404 Cyg has been noted earlier by Jenke et al. (2016), Roques et al. (2015), and Natalucci et al. (2015), who discussed that the high-temperature seed photons might arise from a synchrotron emission in the jet base. However, the radio luminosity during the time of the V404 Cyg preflare observation was very low, ∼10 mJy (Muñoz-Darias et al. 2016; Gandhi et al. 2017), with a steep spectrum that indicates optically thin emission from the jet ejecta (Miller-Jones et al. 2019). The core jet therefore was likely quenched and the synchrotron scenario is not plausible.

The considerations above show why the simple models are not sufficient for modeling the highly unusual hard X-ray spectrum. Natalucci et al. (2015) also fit the hard X-ray data from the INTErnational Gamma-Ray Astrophysics Laboratory with a pure reflection model and a partially covering absorber model. The first model was found to fit most of the data only by assuming very high values for the reflection factor, and the latter by assuming very high column densities for the absorber. The authors regarded these models as unphysical. However, based on the analyses of X-ray and optical data as detailed in Sect. 1, the spectral similarity to the Cyg X–3 hard X-ray state spectrum and the strong X-ray emission lines indicate that the absorption and reprocessing and reflection scenarios are plausible at least for the low-luminosity phases. We therefore continued to fit the V404 Cyg preflare hard X-ray data with absorption and reflection models. We first tried an absorbed power-law model, but found that a fully absorbed model (PHABS × POWERLAW) is completely inadequate to fit the data (${\chi }_{\mathrm{red}}^2=8.2$). Instead, a partially absorbed power law (PCFABS × POWERLAW; model A1 in Fig. 2 and Table 3) fits the hard X-ray data much better (${\chi }_{\mathrm{red}}^2=1.5$), although the intrinsic spectrum is very soft (Γ ∼ 3.9) and partially absorbed (f_cov ∼ 0.93) by a dense medium (N_H ∼ 3.6 × 10²⁴ cm⁻²). An even better fit (${\chi }_{\mathrm{red}}^2=1.1$) can be found by changing the intrinsic spectrum to a cutoff power-law model (A2 in Fig. 2/Table 3) with a value for the power-law photon index more in line with a typical XRB in a hard X-ray state (Γ ∼ 2.2) and cutoff energy of E_cut ∼ 18 keV. The parameters of the absorption component are slightly lower than those of model A1, but they are still comparable. A similar fit quality is achieved by changing the intrinsic emission component to a thermal Comptonization component (model A3 in Fig. 2 and Table 3). The best-fit model has an electron temperature of kT_e ∼ 10 keV, the seed-photon temperature is fixed to 0.1 keV, and the optical depth is τ ∼ 3.6. The absorption component has similar parameter values as models A1 and A2.

For the reflection scenario, we began by fitting the data with PEXRAV (Magdziarz & Zdziarski 1995); a cutoff power-law continuum reprocessed in a neutral medium (model R1 in Fig. 2 and Table 3). We fixed the inclination to the value presented in Table 2 and the abundances to solar. The resulting model has a power-law photon index Γ ∼ 1.8, a cutoff energy E_cut ∼ 20 keV, and a reflection factor R_f ≳ 400. The fit quality is moderate (${\chi }_{\mathrm{red}}^2=1.3$). We also fit the V404 Cyg hard X-ray spectrum with a relativistic reflection model (RELXILL; García et al. 2014; Dauser et al. 2014). The resulting model (R2 in Fig. 2 and Table 3) also has a very high reflection factor (R_f ∼ 90) in a moderately ionized matter (log ξ ∼ 2.3). The intrinsic spectrum is a very hard cutoff power law with Γ ∼ 1.1 and cutoff energy of E_cut ∼ 20 keV. The spin of the black hole is low: a ∼ 0.2. Fixing the inclination to the value presented in Table 2 did not result in a good fit, therefore we left it free to vary, which returned rather low values of θ_inc ≲ 10°. All the other parameters were fixed in the default values. The fit quality is slightly better than in model R1 (${\chi }_{\mathrm{red}}^2=1.1$). In the reflection models, the amount of the radiation that ionizes the reflecting material is typically defined as the ratio of the photon intensity that illuminates it to the direct photon intensity that reaches the observer. When instead of the accretion disk, the reflecting medium is surrounding the X-ray source, the maximum reflection factor can be much higher than in the unobscured case, where usually R_f ≲ 1. We note, however, that high values for the reflection factor, R_f ≲ 10, can be accommodated in the unobscured case as well when strong light bending deep in the gravitational potential of the compact object is assumed (Dauser et al. 2016). Clearly, the reflection factors in the fits are much higher, and we can assume that the majority of the incident emission is reprocessed in a surrounding medium. From the RELXILL model family, we also tried the lamppost geometry (RELXILLLP), which resulted in similar parameters as for model R2 (coronal geometry), but with an even higher reflection factor. Finally, we selected a model with a Comptonized incident spectrum: XILLVERCP (model R3 in Fig. 2/Table 3). Because the reflection factor found above is very high, we fit only the reflection component to the data. The parameters and fit quality are very similar to the cutoff power-law RELXILL model, except for the power-law index, which is Γ ∼ 2. This initial model fitting shows that the absorption or reprocessing models describe the hard X-ray data better than normal continuum models. In addition, when we discard models where the incident continuum is not too soft or too hard for a hard X-ray state, we are left with models A2, A3, R1, and R3.

Next, we included the soft X-ray (3–10 keV) data in the model fitting. Figure 2 (left panel) shows the above models for the whole data range. It is clear that models with partially absorbed but soft intrinsic spectra would need further absorption components to bring the spectrum down to match the data in the soft X-rays, while models with a hard intrinsic spectrum need an additional soft component to account for the data. We selected two models to continue fitting the full data range: an absorbed cutoff power-law continuum (A2), and a fully reprocessed thermal Comptonization continuum (R3). The reasoning behind this selection is that models A2 and A3 are likely very similar, therefore we selected the slightly better fit of model A2. Model R3 was selected because an ionization parameter is included in the model, because of relativistic effects to the spectral shape, and because the fit quality is slightly better.

For the absorbed thermal Comptonization continuum (i.e., model A2), we first added a soft component that we modeled with an absorbed disk blackbody component (PHABS × DISKBB). However, any model producing a Planckian-type spectrum, for instance, thermal Comptonization or bremsstrahlung, produced equally good fits. The absorbed disk with kT ∼ 1.1 keV and N_H ∼ 6.2 × 10²³ cm⁻² can adequately model the soft X-rays, but the soft X-ray slope below 6 keV is not well fit, and large residuals can be found in the 9–10 keV region as well (Fig. 2, middle panel). We therefore added another partial covering absorber and a smeared edge to improve the model. The resulting fit is much better in the soft X-rays, but residuals are still found in the 5–10 keV energy range, likely due to missing line components. Therefore we added two Gaussian lines to the model: an emission line fixed to 6.4 keV (neutral iron), and an absorption line fixed to 6.5 keV (ionized iron). This resulted in a fit quality of ${\chi }_{\mathrm{red}}^2=1.4$. All the parameters of this model can be found in Table 4. A very similar but physically more accurate model was fit to all datasets and is discussed in more detail in the following section.

For the fully reprocessed thermal Comptonization continuum (i.e., model R3), we first fit the same model again for the whole range. The soft X-rays are better taken into account by increasing the value of the ionization parameter from log ξ ∼ 2 to log ξ ∼ 3 (Fig. 2, right panel). Thus, there is no need to add a soft component. Clear residuals are left for the soft X-ray regions below 4 keV and between 7–10 keV. We continued to add an absorption (PHABS) and a smeared edge (SMEDGE) component to bring the model down in these regions. The resulting fit was already much better with ${\chi }_{\mathrm{red}}^2=2.4$, but some small residuals remained below 4 keV and around ∼6.5 keV. Therefore we further added a partial covering absorption component (PCFABS) and an absorption line (GAUSS) to the model, resulting in a fit quality of ${\chi }_{\mathrm{red}}^2=1.4$. A very similar model was fit to all datasets and is discussed in more detail in the following section.

3.3. X-ray spectra: Physical modeling

Because of the observational evidence of a high-density environment described above, we considered the possibility that all sources are embedded in a dense medium (N_H ≳ 10²³⁻²⁵ cm⁻²) that surrounds the X-ray source and causes significant absorption and scattering that affects the X-ray spectrum up to ∼30–40 keV. To facilitate this scenario for spectral fitting, we considered two models: model (A) consisting of a partially absorbed thermal Comptonization component reprocessed in a highly ionized plasma (XILLVERCP, or RELXILLCP when relativistic effects are important) and model (B) consisting of an intrinsic cutoff power-law component (mimicking the thermal Comptonization process) reprocessed in a surrounding neutral uniform-density sphere with polar cutouts of various sizes resembling a torus of different opening angles (BORUS02; Baloković et al. 2018). Model (B) is similar to the Compton-thick AGN scenario where the X-ray source is surrounded by a thick torus with the emission received from a highly absorbed line-of-sight component and a reflected component from the surface of the torus (to follow the discussion of the resulting torus geometries for each source, we refer to Fig. 3). In model (A), the reprocessing occurs in a shell or shells of ionized matter surrounding the X-ray source.

Fig. 3. Geometry of the torus-shaped obscuring matter around the intrinsic X-ray source according to the parameter θtor in model B fits for all sources, in addition to a graphical representation of the column densities (gray, darker for a higher column density), the wind speeds (arrows), and the lines of sight (dashed lines). The thick dashed line around V404 Cyg epoch 1 and 2 corresponds to the partially absorbed blackbody component. The model parameters for each source can be found in Tables 6 and 8.

Model A corresponds to a scenario in which the spectra are dominated by a reflection or scattered component, and the contribution of the incident continuum is severely diminished. This model consists of one or two partially absorbed reflection models (XILLVERCP and/or RELXILLCP), where the incident photons arise from thermal Comptonization (NTHCOMP). Following the indication from the initial modeling in Sect. 3.2 that the majority of incident photons are reprocessed in a medium encompassing the incident photon source, we fixed the reflection factor to a negative value in the model (all radiation was reprocessed). The parameters of the two reflection components (if needed) were kept at the same values, except for the ionization parameter and normalization, which were left free to vary separately for both components. This is to account for changes in ionization parameter in the scattering component, which is evident in the high-resolution X-ray spectra observed from V404 Cyg, Cyg X–3, and GRS 1915+105, which show neutral as well as ionized iron lines. For V404 Cyg and GRS 1915+105, we allowed the redshift to vary freely for the scattering components in order to fit the iron line centroids of ∼6.3 keV, indicating either gravitational or Doppler redshift of the neutral iron line (for Cyg X–3, we fixed this to 1000 km s⁻¹, which is approximately the wind speed of the Wolf-Rayet companion, and to zero for V4641 Sgr because the spectral resolution is far lower). In addition, a narrow iron absorption line and a smeared iron edge with variable absorption energy are needed to successfully fit the first two epochs of V404 Cyg, and epoch 1 of GRS 1915+105. This might indicate an additional absorbing medium during these epochs. Thus, the total model can be described as CONSTANT × PHABS × SMEDGE × PCFABS × (XILLVERCP₁/RELXILLCP + XILLVERCP₂ + GAUSS). Here, constant is the instrument cross-normalization. We did not fix the inclination of the scattering components because the reflecting surface might be inclined away from the disk inclination angle, for example, for an equatorial outflow disk wind with an opening angle of several tens of degrees, or in the case of spherical obscuration, the reflection angle would correspond to some mean angle from all scattering processes. To reduce the parameter space, we fixed the black hole spin to 0. If it is let free, the spectral fits do not restrict the parameter well, which has also been found by Walton et al. (2017) in the case of V404 Cyg, where the spin was estimated as a >  −0.1. It can be also expected that the reflecting medium can lie farther away than the innermost stable circular orbit (ISCO) when reprocessing in the surrounding media is assumed. We therefore fixed the outer radius to 1000 gravitational radii in all sources.

For model B, we used the BORUS02 model component (Baloković et al. 2018), which instead of a disk reprocessing allows a variety of geometries from a uniform sphere to torus-like shapes through polar cutouts. The reprocessing torus in BORUS02 is considered to be cold, neutral, and static. To take a moving reprocessor into account, we therefore allowed the redshift of the scattering component to vary. Like in model A, we also introduced a highly ionized iron edge component and an ionized iron absorption or emission line component to find acceptable fits to the data. We can expect that the reprocessed photons are either redshifted (behind the source, i.e., moving away) or blueshifted (in front of the source, i.e., moving toward) because they arise in the fast outflow or wind or in the accretion flow. The total model can be described as CONSTANT1 × PHABS1 × SMEDGE1 × (CONSTANT2 × BORUS02(RED) + CONSTANT3 × BORUS02(BLUE) + PHABS2 × CABS1 × CUTOFFPL + CONSTANT4 × CUTOFFPL). CONSTANT1 is the instrument cross-normalization, CONSTANT2 is the relative normalization of the redshifted scattered component, CONSTANT3 is the relative normalization of the blueshifted scattered component, CONSTANT4 is the relative normalization of the leaked (unabsorbed) intrinsic spectrum, PHABS2 × CABS1 is the line-of-sight absorption including beam-scattering with the column densities linked between the components, and CUTOFFPL represents the intrinsic continuum of the accretion flow (mimicking thermal Comptonization spectrum). For the scattered component, we fixed the inclination angle according to Table 2 because the scattering angle is now taken into account in the model. We also fixed the iron abundance to solar. All the parameters are linked between BORUS02(RED) and BORUS02(BLUE), except for the redshift for BORUS02(BLUE), which is determined as being the negative of the value in BORUS02(RED). In addition, for the first two epochs of V404 Cyg, an additional soft component is needed, which we modeled with a partially absorbed blackbody component. However, any component resembling a Planckian spectrum might be inserted instead, such as a low-temperature thermal Comptonization component (cf. model A), a disk blackbody component, or a Wien spectrum. The physical interpretation of this component is discussed in Sects. 3.3.2 and 4.1. In addition, as mentioned above, we added an emission line component to the preflare spectrum of GRS 1915+105 and to the Cyg X–3 spectrum. In the following, we concentrate on the fitting results of these two models for individual sources.

3.3.1. Cyg X–3

The peculiarity of the hard-state spectrum in Cyg X–3 has been known for more than a decade (Hjalmarsdotter et al. 2004). The observed low-energy cutoff has previously been attributed to either strong absorption in the stellar wind, pure Compton reflection in a medium that covers the emitting source, or nonthermal Comptonization of a steep electron population by Hjalmarsdotter et al. (2008). They preferred the latter scenario, although the former two produce better fits and require either an unusual accretion state or a very massive black hole as a primary star. Later, Zdziarski et al. (2010) showed that a low-energy cutoff in the X-ray spectrum can be obtained when Compton downscattering is considered in an optically thick plasma cloud, likely arising from the interaction of the accretion disk and the strong stellar wind of the Wolf-Rayet companion. In addition, there is evidence that the primary star in Cyg X–3 is a black hole with a relatively low mass (Zdziarski et al. 2013; Koljonen & Maccarone 2017). Because of the different types of companion stars in V404 Cyg, V4641 Sgr, GRS 1915+105, and Cyg X–3, it seems unlikely that the similar X-ray spectra would be a result of different accretion mechanisms; wind versus Roche-lobe accretion. In addition, the hard X-ray state of Cyg X–3 is a relatively stable state that can last up to several years. Therefore it seems unlikely that it would present a peculiar accretion state for such a long time. Rather, the spectral similarity to V404 Cyg, V4641 Sgr, and GRS 1915+105 likely arises from some type of radiation reprocessing.

The interstellar absorption for Cyg X–3 is relatively high as a result of its location in the plane of the Galaxy, and likely because it is located behind two spiral arms and the Cygnus X star-forming region (McCollough et al. 2016). We fixed the lower limit of the hydrogen column to 3.5 × 10²² atoms cm⁻², which is approximately the value found in studies where instruments with softer X-ray response were used (Koljonen et al. 2018; Kallman et al. 2019). Because the reprocessing matter in the model is neutral, we added a smeared edge in model B to account for the absorption of highly ionized iron. The energy of the smeared-edge component is about 9 keV, indicating that the Fe XXVI Lyα line (with an ionization energy of 9.2 keV) is the most dominant absorber. We found that including a narrow line at the energy of 7.8 keV (either from Fe XXV Kβ and/or Ni XXVII Kα) improves the fit as well. We also fixed the redshift in both models to 1000 km s⁻¹ (z = 0.003), which is approximately the wind speed of the Wolf-Rayet companion (Koljonen & Maccarone 2017), although this has only a slight effect on the fit.

The resulting parameters for model A fits can be found in Table 5. The Cyg X–3 hard-state spectrum can be adequately described by a rather soft, thermally Comptonized (Γ ∼ 2.4, kT_e ∼ 33 keV) spectrum reprocessed in a highly ionized medium (log ξ ∼ 3.8) and further in a lower ionization medium (log ξ ∼ 2.9). This probably is a scattering cloud of decreasing ionization parameter. Similar modeling with two ionized absorbers was successfully used in Kallman et al. (2019) to fit the X-ray emission lines from Chandra data with a medium-ionization component having log ξ ∼ 2.9 and a high-ionization component fixed to log ξ ∼ 4.2−5.0. The emission is further partially absorbed in a relatively dense environment (N_H ∼ 4 × 10²³ atoms cm⁻²) with a covering fraction of f_cov ∼ 0.5. This is consistent with the estimates of the wind column of ∼10²³ atoms cm⁻² acquired from fitting the X-ray emission lines (Kallman et al. 2019). Furthemore, the covering fraction can be understood as half of the emission back-scattering and passing through the cloud.

Model B (parameters can be found in Table 6) delivers similar results, with the surrounding dense torus essentially a sphere (cos(θ_tor) > 0.96, N_H, tor ∼ 10²⁴ atoms cm⁻²; see also Fig. 3). The incident spectrum is very similar to that of model A (Γ ∼ 2.4), with a cutoff energy at 28 keV (although if the mechanism is thermal Comptonization, this corresponds to an electron temperature of kT_e ∼ 9−14 keV). The line-of-sight component is absorbed by a column of N_H ∼ 3 × 10²³ atoms cm⁻² and comprises 65% of the flux received, while the 17% and 18% of the flux come from the scattered and leaked (or direct) component, respectively. The unabsorbed 3–79 keV luminosity is similar for both models and corresponds to ∼10% of the Eddington luminosity¹ for a 2.5 solar mass black hole, which is fairly high and corresponds to values found for V404 Cyg. However, because Cyg X–3 is orbiting the Wolf-Rayet star inside its photosphere, there is always ample matter to accrete and sustain the high luminosity. On the other hand, assuming a 10 solar mass black hole (an upper limit for the allowed mass of the compact object; Koljonen & Maccarone 2017), the luminosity would be ∼2–3% of the Eddington luminosity, which is similar to what has been observed from Cyg X–1 (Basak et al. 2017).

3.3.2. V404 Cyg

V404 Cyg (GS 2023+338) is one of the closest XRB (2.39 ± 0.14 kpc; Miller-Jones et al. 2009) with a relatively long orbit (6.473 ± 0.001 d; Casares et al. 1992), implying a large accretion disk and a long outburst recurrence time. Overall, two outbursts have been detected from V404 Cyg with X-ray instruments in 1989 (Kitamoto et al. 1989; Oosterbroek et al. 1997; Życki et al. 1999) and 2015 (Rodriguez et al. 2015; Motta et al. 2017a, b; Sánchez-Fernández et al. 2017; Walton et al. 2017; Kajava et al. 2018; considering the June and December 2015 flaring events to be parts of the same outburst), although in retrospect, additional optical outbursts have been detected in 1938 and 1956 (Richter 1989; Wagner et al. 1991). The two X-ray outbursts in 1989 and 2015 were both hard-state outbursts with no excursion into a soft X-ray state, although the luminosities reached or exceeded the Eddington luminosity. During the 2015 outbursts, a plethora of X-ray spectral behavior was observed from V404 Cyg, with some applicable to absorption events with intermediate flux densities and X-ray spectra not consistent with a Comptonization model, and some to intrinsic variations with very high or low flux densities consistent with a Comptonization model (Motta et al. 2017a, b; Sánchez-Fernández et al. 2017; Kajava et al. 2018; Walton et al. 2017; Hynes et al. 2019). In both cases, the X-ray spectra exhibit fast changes from one state to another in a matter of seconds to minutes (Motta et al. 2017a; Kajava et al. 2018; Sánchez-Fernández et al. 2017; Walton et al. 2017).

We identify three different spectral states in the set of 30 000 cts spectra during the NuSTAR pointing (Figs. 4 and 5). At the beginning of the pointing, the source presents Cyg X–3 hard-state-like spectra, that is, spectra that likely exhibit strong absorption or reprocessing, as evidenced by the curved hard X-ray spectra (epoch 1 in Fig. 4, blue spectra in Fig. 5). This spectrum was used in the initial model fitting in Sect. 3.2. The spectrum further evolves to a softer state with an increase in the fluxes below 6 keV and a decrease in the fluxes above 10 keV (epoch 2 in Fig. 4, yellow spectra in Fig. 5). Both epochs are approximately 10 ks long, and the change takes place in a gap between spectra 6 and 7 (Fig. 6). During the transition, the energy region between 7 and 10 keV remains remarkably constant. This might be interpreted as the soft and hard part of the spectrum arising from two different components that are anticorrelated (spectral pivoting is not enough to explain the whole spectral change). After epochs 1 and 2, V404 Cyg entered into a high count rate flaring period and presented significantly harder spectra and a variable cutoff energy higher than that of the preflare spectra (epochs 3–5 in Fig. 4, green spectra in Fig. 5).

Fig. 4. NuSTAR 3–79 keV light curve of V404 Cyg divided into GTIs that total 30 000 counts (numbered). The numbers for GTIs that include count rates exceeding 1050 cts s−1 are not shown (excluding GTIs 13 and 14), and the data points are shown in gray with a reduced point size. The coloring scheme corresponds to the spectral shape shown in Fig. 5, with every other GTI shown in a different hue for clarity.

Fig. 5. Changes in the spectral shape of V404 Cyg during the GTIs shown in Fig. 4. All the spectra have been normalized to the 7 keV flux in the first spectrum to show the difference in the spectral shape.

Fig. 6. Spectra of V404 Cyg from GTIs 6–9 in Fig. 4, showing the change from a harder (6) to a softer spectrum (7–9) preceding the flaring period. The spectra are not normalized, thus showing the variable regimes below 6 keV and above 10 keV.

We fit the averaged spectra from epochs 1–5 and the individual GTIs with models A and B. The model parameters, corresponding fluxes, and the fit quality for the averaged spectra can be found in Table 7 for model A and in Table 8 for model B, while a selection of the model parameters and fluxes is shown for the individual GTIs in Figs. 7 and 8 (with a fit quality ranging between ${\chi }_{\mathrm{red}}^2$ = 0.9–1.6 with a mean of ${\chi }_{\mathrm{red}}^2$ = 1.2). In addition, the average spectra and the corresponding model B fits divided into different model components (both absorbed and unabsorbed) are shown in Fig. 9.

Fig. 7. Model A parameter values for V404 Cyg spectra (individual GTIs in black, epochs in color).

Fig. 8. Model B parameter values for V404 Cyg spectra (individual GTIs in black, epochs in color).

Fig. 9. Averaged NuSTAR FPMA data from V404 Cyg in the preflare (epochs 1–2) and outburst (epochs 3–5) stages as shown in Fig. 5, together with absorbed (top row) and intrinsic (bottom row) model B components. The total model (solid red line) consists of a sum of a blackbody component (dashed green line), a cutoff power-law component absorbed in and scattered off from the material in the line of sight (dotted blue line), scattered into the line of sight by the surrounding medium (dot-dashed light blue line), and unabsorbed or direct component (long-dashed yellow line). The middle panels show the residuals of the models to the data.

Considering model A fits, the spectra of all epochs can be fit with the fully reprocessed thermal Comptonization emission; a model that is essentially very similar to the basic model used in Walton et al. (2017, their Table 3), except that here we used the version of RELXILL with a spherically symmetric corona and an intrinsic emission arising from thermal Comptonization (RELXILLCP) instead of a lamp-post geometry and an intrinsic emission modeled as a cutoff power-law spectrum (RELXILLLP). In Sect. 3.2 the lamppost model fit to the epoch 1 hard X-ray spectrum (essentially very similar to the fit with coronal geometry) would require a very steep incident power-law spectrum. In addition, we used a simple absorption component instead of an XSTAR reprocessor. We also let the ionization parameter and the redshift of the XILLVERCP component vary freely instead of fixing them to 1 and 0, respectively, to fit the redshifted iron line with energies ∼6.3 keV. The results of the model A fits are very similar to those in Walton et al. (2017) for epochs 3–5 (the flaring state), as expected, with low intrinsic absorption, an inclination close to 30°, a photon index of ∼1.6 on average, and an ionization parameter of ∼1000. Some differences do arise, however, with the iron abundance close to solar in our fits except for a few GTIs in epoch 4, compared to twice the solar value in Walton et al. (2017). On another note, our fits for the epoch 1–2 averaged spectra also show elevated abundances (for individual GTIs, the abundances were fixed to the averaged value of the epoch because they were not well constrained in the fits). Because the geometry and intrinsic emission were modeled differently, the remaining parameters are more difficult to compare, but Walton et al. (2017) reported very high values for the reflection factor (R_f ∼ 1−3). This might also indicate reprocessing in the surrounding medium, as speculated in this paper, instead of a strong gravitational bending and a scenario with a high black hole spin.

On the other hand, epochs 1–2 display much softer spectra, as shown in Fig. 5, which is reflected in the fits with increased values for the photon index (Γ ∼ 2) and ionization parameter (log ξ ∼ 3.5), and very low values for the electron temperature (kT_e ∼ 3−4 keV, corresponding to an optical depth of τ ∼ 10). In addition, a partially covered absorption component with a high column density (4 × 10²³ cm⁻², f_cov ∼ 0.8), a narrow iron absorption line at 6.5 keV, and a smeared iron edge at 7.4 keV likely arising from an additional variable absorption component are needed to fit the spectra successfully with RELXILLCP (an additional XILLVERCP component is not needed in these epochs). The physical explanation for the low electron temperature is a challenge in this model. Because the cutoff energy in the epoch 1–2 spectra is low (< 20 keV), all models with thermal Comptonization are expected to give electron temperatures that are approximately lower than 6–10 keV (assuming E_cut ∼ 2−3 kT_e). We recall that the electron temperature in the model is given in the frame of the observer. A very efficient cooling mechanism, such as radiative cooling by soft photons from the strong outburst or reprocessing in an optically thick medium, can therefore thermalize and Compton downscatter the intrinsic emission to lower energies (similar to what has been proposed for Cyg X–3 in Zdziarski et al. 2010, see also Sect. 3.3.1). The increase in soft X-ray emission between epochs 1 and 2, as mentioned above, is mirrored in the model A parameter evolution as a change in the absorption parameters (decrease in the column density and covering fraction), as an increase in the power-law photon index and normalization, and a decrease in the electron temperature. In addition, there is an increase in R_in, log ξ, z, and inclination. In Sect. 4.1 we discuss the hypothesis that this parameter evolution can arise from a geometry change gearing toward jet ejection.

In epochs 1–2, model B fits are consistent with the torus completely covering the source (cos(θ_tor) pegged to 1; see also Fig. 3), with an average density of the torus of N_H, tor ∼ 10²⁴ atoms cm⁻² and a similar if slightly higher (a factor of two) line-of-sight column density through the torus. The flux from the line-of-sight component dominates the spectrum during these epochs, comprising of 75–90% of the total flux, while scattered and leaked flux contribute only up to 2% (the remaining flux comes from the blackbody component discussed below). In epochs 3–5, the geometry changes from spherical to more disk-like with cos(θ_tor)∼0.55, and an increase of a factor of four in N_H, tor, while the line-of-sight column density decreases by a factor of about six (see also Fig. 3). The flux from the line-of-sight component decreases to 25–50% of the total flux, while the scattered flux increases to 40–70% of the total flux and causes the majority of the spectral variability in this state. Moreover, the flux from the leaked component increases significantly, up to 10% of the total flux, indicating that the geometry has changed so that the direct intrinsic emission is partly visible. Similar to the parameter evolution in model A, the intrinsic cutoff power-law spectrum is soft in epochs 1–2, with the power-law photon index Γ pegged to 2.6 and cutoff energy values of about 21–26 keV. The spectrum hardens in epochs 3–5 with Γ ∼ 2 and much higher cutoff energies, although overall the intrinsic spectrum is softer than in model A. In contrast to model A, no absorption line is required in the model for epochs 1–2 because the double-horned line profile can be explained by reprocessing in an outflowing medium, which is seen as both redshifted and blueshifted (z ∼ 0.05c ∼ 15 000 km s⁻¹). Previously, Motta et al. (2017a) have shown that the iron line exhibit redshift and blueshift with velocities v <  0.1c, in line with the NuSTAR values. For epochs 3–5, the resulting line speeds (∼4000 km s⁻¹) are comparable to the values measured from the P Cygni profiles of He I λ5876 (a well-known accretion wind tracer), while for epochs 1–2, they are still well within of what can be achieved with super-Eddington accretion (usually on the order of 0.1c − 0.2c; e.g., Pinto et al. 2020). In epochs 3–5, the back-illumination dominates the scattering (i.e., c2 > c3 in Table 8), so that the scattered emission is seen through material that is more translucent than that of the line of sight. The scattered line therefore exhibits mainly redshift and reduces the observed iron Kα line energy to 6.3 keV without producing the blueshifted line.

The strong absorption of the intrinsic emission in model B requires an additional component to account for the elevated soft X-ray emission in epochs 1–2. We modeled this as a partially absorbed blackbody component, as discussed in Sect. 3.3. The temperature of the blackbody component decreases from 1.1 keV to 0.95 keV, while its normalization (and flux) increases from epoch 1 to epoch 2, indicating an increase in the size of the emitting medium. While the temperature of this component might arise from the hot disk, it does not seem plausible that the emitting area would increase and temperature decrease when the source is gearing toward the flaring state when emission from an accretion event is considered. Życki et al. (1999) briefly discussed that the blackbody component needed to provide the soft excess in the 1989 outburst might arise from incident disk photons that are thermalized in multiple scattering processes in the surrounding medium. Thus, a dense medium, whether a disk with a large scale-height or a stellar or accretion disk wind, might then explain both the large column and the soft component, similar to what was discussed for model A above.

For both models, the unabsorbed luminosities in the 3–79 keV NuSTAR band are sub-Eddington for all epochs (2–4% for model A, and 2–14% for model B), although because the intrinsic spectrum is steep in epochs 1–2, the Eddington limit would be reached by extrapolating the spectrum down to ∼0.2 keV (for model B at least). The luminosity of epochs 1–2 is much higher for model B because the highly absorbed intrinsic emission is included, while model A presents only the flux for the scattered component. When the scattered flux is only taken into account in model B, the resulting luminosities agree with those for model A.

3.3.3. V4641 Sgr

V4641 Sgr is a very peculiar source, presenting a dynamically confirmed black hole with a high-mass companion (Orosz et al. 2001), but showing transient outbursts similar to low-mass XRBs. The outbursts of V4641 Sgr can be very short (about a week) and intense, reaching super-Eddington levels, as in 1999 outburst (Revnivtsev et al. 2002), or longer but much weaker (e.g., Uemura et al. 2002; Maitra & Bailyn 2006; Muñoz-Darias et al. 2018). There is some evidence that the inner accretion disk is misaligned to the orbital plane of the binary (i ∼ 70°; Orosz et al. 2001; MacDonald et al. 2014; Pahari et al. 2015), and V4641 Sgr might instead be a low-inclination source in X-rays and radio (i ∼ 10°; Hjellming et al. 2000; Orosz et al. 2001; Gallo et al. 2014). It has previously been suggested that a Chandra spectrum observed during an outburst decline is very similar to those observed from Seyfert-2 AGN (Morningstar et al. 2014).

The RXTE/PCA spectra of V4641 Sgr gathered from the 1999, 2002, 2003, and 2005 outbursts are plotted in Fig. 10. In the top panel, the spectra from 1999 and the latter part of 2003 are very reminiscent of the spectra with low absorption and high count rate of V404 Cyg that were studied in Walton et al. (2017), while in the bottom panel, the spectra from the 2002, 2003, and 2005 outburst are more reminiscent of the flaring spectra with lower count rate of V404 Cyg (Fig. 5). The only outburst containing both types of spectra is the 2003 outburst, which began with those presented in the bottom panel (epochs 1 and 2 in Maitra & Bailyn 2006, including the spectra shown in Fig. 1) and continued with those presented in the top panel (epochs 3 and 4 in Maitra & Bailyn 2006).

Fig. 10. Changes in the RXTE/PCA spectral shape of V4641 Sgr during the outburst periods of 1999, 2002, 2003, and 2005. All the spectra have been normalized to the 7 keV flux in the first spectrum to show the difference in spectral shape. Top panel: outburst spectra reminiscent of the spectra with high count rate and low absorption of V404 Cyg (Walton et al. 2017), bottom panel: spectra reminiscent of the flaring spectra with lower count rate and higher absorption of V404 Cyg in Fig. 5.

The two models (A and B) were fit to the spectra from two pointings observed during the 2002 (hereafter epoch 1) and 2003 (hereafter epoch 2) outbursts (Fig. 1; the top corresponds to epoch 1, and the bottom corresponds to epoch 2), and the resulting parameters are shown in Tables 5 and 6. Both observations were taken during the middle of the outburst peak with similar optical magnitudes (∼11.5 mag in V and R band; Uemura et al. 2002; Maitra & Bailyn 2006) and X-ray count rates. Because the statistics of RXTE/PCA are much lower than those of NuSTAR, we fixed the inclination to i = 70° in both models to reduce the number of free parameters. In the case of model A, a successful fit could be achieved with a single RELXILLCP model, while for model B, the spectra are dominated by the (back-)scattered component (75–80% of the total flux) for both epochs. The main difference in the parameters of both models for the two epochs can be found in the intrinsic spectrum with the power-law photon index being higher for epoch 1 (Γ_A = 2.2, Γ_B = 2.6) with no high-energy cutoff, while for epoch 2, the power-law photon index is much lower (Γ_A = 1.5, Γ_B = 1.6) and the spectrum has a low-energy cutoff (kT_e ∼ 9 keV, E_cut ∼ 20 keV).

Unlike in other sources, the covering fraction and the line-of-sight absorption column in model B are low for both epochs. The scattered component dominates the flux in model B, while pure scattering in a plasma with a single-ionization parameter is consistent according to model A. In contrast to V404 Cyg epochs 1–2 and GRS 1915+105 epoch 1 (discussed in the next section), there is evidence of a strong jet in the two V4641 Sgr epochs. The 8.5 GHz radio observations of the Very Large Array show flux densities of 80–170 mJy (Rupen et al. 2002) and 550–570 mJy (Rupen et al. 2003) coinciding with epochs 1 and 2, respectively. This means that during both observations, any surrounding matter was likely evacuated by the jet. Because the inclination of the system is likely high, we might observe the X-ray source through a disk wind or geometrically thick accretion flow (the column density of the torus is the highest of all the four sources for model B fits; log N_H, tor ∼ 10²⁵ cm⁻²), and most of the emission received is from backscattering (i.e., c2 = 0 in Table 6; see also Fig. 3).

The luminosities of the two V4641 Sgr epochs are much lower (∼5 × 10³⁶ erg s⁻¹, corresponding to ∼0.5% of the Eddington luminosity) than the luminosities for the other sources, although the observations are from the peak of the outbursts. The 2002 and 2003 outbursts were weaker than the much more luminous outburst in 1999, where the peak luminosity was at or greater than the Eddington luminosity (Hjellming et al. 2000; Revnivtsev et al. 2002). These weak outbursts do not have any predictable periodicity, but they seem to occur roughly at intervals of 500–600 days (Negoro et al. 2018). On the other hand, an optical counterpart as bright as in the 2002 and 2003 outbursts has not been detected for any other weak outbursts since the 1999 outburst, which marks them as different and probably means that they included enhanced reprocessing of the X-ray emission in the accretion disk.

3.3.4. GRS 1915+105

GRS 1915+105 is one of the brightest XRBs in our Galaxy in its outburst because it has the longest orbital period known among low-mass XRBs (33.9 days; Steeghs et al. 2013) and thus the largest accretion disk size because of the largest tidal truncation radius. The huge mass reservoir has lasted already three decades, powering the outburst until the drop in the X-ray flux in 2018. We fit all three NuSTAR observations (labeled epochs 1–3) from this anomalous state with models A and B. The resulting model parameters are shown in Tables 5 and 6. With model A, epoch 1 is fit best with a single RELXILLCP model, while epochs 2 and 3 can be fit with two RELXILLCP models, one corresponding to a plasma with a higher (log ξ ∼ 3.4) and the other to a plasma with a lower ionization parameter (log ξ ∼ 2). Epoch 1 includes a prominent narrow (up to detector resolution) absorption line at 6.56 keV that might be the Fe XXV Kα line, although it is redshifted by 6000–7000 km s⁻¹, or the Fe Kα line blueshifted by the same amount. In addition, an edge at 7.4 keV is needed to fit the spectra adequately. In the case of blueshifted neutral iron absorption, this implies a velocity of the absorbing material of ∼13 000 km s⁻¹. However, the redshift needed to model the spectrum is even higher, z ∼ 0.12, which agrees with what is observed from ultraluminous X-ray sources (ULXs). For epochs 2–3, the redshift required is much lower (∼3000 km s⁻¹). In contrast, the iron abundance increases from solar in epoch 1 to around twice the solar value in later epochs. The incident cutoff power-law spectrum changes from epoch 1 to epochs 2–3, displaying a decrease in the power-law photon index from 1.9 to 1.6, with the latter two also showing a cutoff at 14 keV, while in epoch 1, the cutoff is unconstrained (≳160 keV, when left free to vary). The inclination is only constrained for epoch 1, corresponding to 60–80 degrees, while for epochs 2 and 3, we fixed it to 70°. The models are absorbed with N_H ∼ 5×10²² atoms cm⁻², which is comparable with the interstellar value of N_H ∼ 3.5 × 10²² atoms cm⁻². Additional, partial absorption is needed for epochs 2 and 3, with N_H ∼ 5 × 10²³ atoms cm⁻² and a covering fraction of ∼0.5.

With model B, epoch 1 is consistent with the torus covering the source completely, while cos(θ_tor) decreases to 0.8 and 0.7 in epochs 2 and 3, respectively (Fig. 3). The density of the torus remains equal in all epochs (N_H, tor ∼ 2.5 × 10²⁴ atoms cm⁻²). The line-of-sight component contributes about equally in all epochs (52–59%) with a similar column density (N_H, los ∼ 2.5−5 × 10²³ atoms cm⁻²). The scattered flux increases from 17% to 32–35%, while the leaked flux decreases from 31% to 9–10% when transiting from epoch 1 to epochs 2 and 3. The scattering component is dominated by the redshifted scattering, with epoch 1 showing high wind speeds of 0.057c (∼17 000 km s⁻¹, which is more in line with the velocity of the absorption line if it is due to a blueshifted neutral iron line) that reduce to much lower values (600–2700 km s⁻¹) in epochs 2 and 3. The incident spectrum in epoch 1 is a steep power law (Γ ∼ 2.3) similar to V404 Cyg epochs 1, 2, and 5, and Cyg X–3 spectra. In epochs 2 and 3, the power-law index of the incident spectrum decreases to 1.6–1.8 and the spectrum exhibits a cutoff at 21–24 keV. Interestingly, the unabsorbed flux remains similar in all epochs and corresponds to 1–2% of the Eddington flux. The average spectra and the corresponding model B fits divided into different model components (both absorbed and unabsorbed) are shown in Fig. 11. The inclination is fairly well constrained in epochs 2 and 3 and corresponds to 40–60°. For epoch 1, it is not well-constrained in the fit, and we froze it to 53° (cos(θ_inc) = 0.6). Based on measurements of the jet inclination of the system, Reid et al. (2014) estimated the disk inclination angle as 60° ± 5°, which is consistent with the values derived here.

Fig. 11. Averaged NuSTAR FPMA data from GRS 1915+105 from the three epochs together with absorbed (top row) and intrinsic (bottom row) model B components. The total model (solid red line) consists of the sum of a cutoff power-law component absorbed in and scattered off from the material in the line of sight (dotted blue line), scattered into the line of sight by the surrounding medium (dot-dashed light blue line), and unabsorbed or direct component (long-dashed yellow line). Middle panels: residuals of the models to the data.

While having slightly different parameters, the two models show a similar evolution of the model parameters, consistent with a scenario of obscured emission through fast (spherical) outflowing wind in epoch 1, which flattens and decelerates in later epochs (Fig. 3). Epoch 1, which is observed immediately after the sudden decrease of the hard X-ray emission, shows higher wind velocities but neutral absorbing material, higher incident power-law photon indices, and no spectral cutoffs. This is consistent with a scenario where the inner accretion flow is obscured by matter that is first seen as scatterer and absorber. Later, it is mostly seen in reflection, and the regions with higher ionization are exposed.

3.4. X-ray timing

While we did not concentrate on the X-ray timing properties of the sources in detail, in the following we provide a quick analysis of the NuSTAR X-ray cospectra (a proxy for the PSD). The cospectra for the NuSTAR data for all sources are shown in Fig. 12. As described in Sect. 2, we extracted two cospectra for the soft X-ray band (3–10 keV) and the hard X-ray band (10–79 keV), and when they did not present any differences from each other, we extracted the cospectra from the full NuSTAR range 3–79 keV.

Fig. 12. Cospectra of V404 Cyg, Cyg X–3, and GRS 1915+105 for NuSTAR epochs. The soft- and hard-band cospectra are shown separately when they display any differences. The best-fit single or broken power-law model is also plotted for all data. Additional Lorenzian functions are needed for the QPO and harmonics in the GRS 1915+105 epoch 1 data.

We can assume that the reprocessing or scattering in the obscuring matter smears out most of the high-frequency timing information of the intrinsic emission. For Cyg X–3, it has been shown that the PSD is close to a power law with an index of −2.0, independent of the accretion state of the source (Axelsson et al. 2009; Koljonen et al. 2011), and this has been speculated to be due to a suppression of the high-frequency variations by scattering in the stellar wind surrounding the X-ray source, mimicking a red noise process (Koljonen et al. 2018). Likewise, for the NuSTAR pointing considered here, we obtain a power-law cospectrum with an index of −1.9 ± 0.1.

For V404 Cyg, a timing analysis has been performed on a subset of the NuSTAR data used in this paper in Gandhi et al. (2017). The authors found that the X-ray cospectrum is consistent with a power-law spectrum with an index of −1.6; this is not as steep as in Cyg X–3 or V4641 Sgr, but steeper than the flicker noise that is typically observed from hard states of XRBs, indicating the suppression of high-frequency variability. We studied here the epoch-by-epoch variations in the cospectrum and found that a single or broken power law fits all the data sufficiently well. In all epochs, the slope of the high-frequency spectra agrees roughly with an index of ∼−1.6 (epoch 1: −2.0${}_{-0.2}^{+0.4}$/−1.41 ± 0.06, epoch 2: −1.58 ± 0.08/−1.41${}_{-0.03}^{+0.08}$, epoch 3: −1.6${}_{-0.1}^{+0.2}$, epoch 4: −1.7 ± 0.1, epoch 5: −1.4 ± 0.1; when two numbers are given, they correspond to the 3−10 and 10−79 keV band), while there is some evidence of a break to flatter indices at lower frequencies. Interestingly, the epoch 2 cospectrum shows a diminished rms in the soft X-ray band compared to the hard X-ray band, which might incidate that the soft spectral component dilutes the rms. In addition, the spectral break shifts to lower frequencies from 0.2 Hz to 0.06 Hz, and if this is related to the size of the varying soft component, it indicates an increasing emitting area that is consistent with the blackbody parameter evolution of model B fits.

GRS 1915+105 is famous for its complex X-ray variability. The PSDs typically consist of a band-limited noise component with one or more peaks, indicating quasi-periodic oscillations (QPOs; e.g. Morgan et al. 1997). The cospectra from epochs 1–3 are consistent with a power-law noise (epoch 1: 1.21 ± 0.02/0.91 ± 0.05, epoch 2: 1.02 ± 0.06, and epoch 3: 1.96 ± 0.08; when two numbers are given, these correspond to the 3–10 and 10–79 keV band), with the epoch 1 cospectrum also showing a low-frequency QPO with harmonics at 3.3 Hz (at least for the soft X-ray band). Clearly, the epoch 1 cospectrum differs from epochs 2–3, which do not present power over 0.1 Hz, and likely we have a more direct view of the accretion flow. In addition, the soft X-ray band where the incident (leaked) spectral component dominates (Fig. 11) is more variable than the hard X-ray band. For the model B fit, the leaked component presents the highest luminosity fraction in epoch 1, comprising 31% of the total flux, and in this pointing, we might be observing a patchy outflow. For epochs 2–3, the amount of flux observed from the leaked component is diminished and the total flux is dominated by the line-of-sight and scattered components, which accounts for the loss of high-frequency power. The low-frequency power varies from flicker noise in epochs 1–2 to red noise in epoch 3.

For V4641 Sgr, Maitra & Bailyn (2006) studied the X-ray timing properties during the 2003 outburst and found red-noise-dominated PSD below 1 Hz and Poisson noise above it. Slightly more structure was found in the much brighter 1999 outburst, with the PSD showing a broken power-law shape with indices of ∼−1 and ∼−2 below and above 5 Hz, respectively (Wijnands & van der Klis 2000). However, no other structure, such as QPOs, were found in the PSD reaching 100 Hz, and again the higher frequency spectrum is consistent with the red-noise process. For the RXTE pointings of V4641 Sgr considered here, the PSDs are consistent with pure Poisson noise in the case of the 2002 data, and a red noise below 0.03 Hz and Poisson noise above for 2003 data.

4. Discussion

We have shown that the X-ray spectra of V404 Cyg, Cyg X–3, GRS 1915+105, and V4641 Sgr share similarities at certain evolutionary times that correspond to low X-ray flux periods during outburst events (considering that the persistent source Cyg X–3 is always “on”, that is, in outburst, and that V4641 Sgr exhibits low-luminosity outbursts). We modeled these spectra successfully with two models consisting either of a fully reprocessed spectral component and/or a heavily absorbed spectral component, with the intrinsic spectra arising from a thermal Comptonization process. It is well known that Cyg X–3 orbits its Wolf-Rayet companion star inside a high-density stellar wind and that it can strongly affect the X-ray spectra. Based on the spectral similarity, we can therefore presume that similar surroundings affect the X-ray spectra of V404 Cyg and V4641 Sgr during the outbursts, and the recent anomalous accretion state of GRS 1915+105. Because the companion stars of V404 Cyg, V4641 Sgr, and GRS 1915+105 are not expected to present high-density stellar wind, we attribute this medium to either a large scale-height accretion flow (possibly due to a super-Eddington accretion rate) or an optically thick equatorial outflow or envelope either from the radiation pressure of intrinsic (super-Eddington) emission of the accretion flow or from the base of the jet.

Both models (A and B) indicate a change in geometry in the system in the evolution of V404 Cyg and GRS 1915+105. Epochs 1–2 for V404 Cyg and epoch 1 for GRS 1915+105 are consistent with the X-ray source being obscured in an outflowing (spherical) plasma cloud that transforms into a more disk-like geometry in subsequent epochs (Fig. 3). This change also coincides with the start of the activity in the radio and X-ray flaring. The two V4641 Sgr epochs are consistent with a disk-like geometry, and in both cases, strong radio emission was detected at the same time. This means that the geometry change of the obscuring component seems to be linked with the radio evolution. Interestingly, a very similar evolution took place in the 1999 outburst of V4641 Sgr (Revnivtsev et al. 2002), with the source luminosity dropping by an order of magnitude followed by optical and subsequently radio emission and a change in the X-ray spectrum from softer (epoch 1-type) to harder (epoch 2-type).

In V404 Cyg and GRS 1915+105, the intrinsic spectra first resemble that of an intermediate state of XRBs with power-law photon indices Γ ∼ 2.0−2.6, and changes to lower values, Γ ∼ 1.4−2.0, more typical of an XRB hard-state spectrum. However, their accretion history is very different. V404 Cyg epochs 1–2 are the softest state that the source goes through during the outburst, while for GRS 1915+105, epochs 1–3 are the hardest so far observed from the outburst (epoch 3 being much harder than epochs 1–2). Therefore it does not seem likely that a single accretion state would explain the similar-looking spectra. Rather, it has to do with the reprocessing of the intrinsic emission. In the following, we discuss this issue further.

4.1. Was V404 Cyg in a “soft-state” during epochs 1–2?

The outburst spectra of V404 Cyg do not show a soft blackbody emission component, except for the 1989 outburst, when a 0.3 keV disk component was seen; see Życki et al. (1999). We discussed in Sect. 3.3.2 that the parameter evolution of the blackbody component is not consistent with arising directly from the accretion disk. Instead, the soft component could arise from incident thermal (disk) photons scattered multiple times in the surrounding changing medium. On the other hand, the fits with the reflection model in Sect. 3.3.2 showed that the soft component could also be modeled with an increasing ionization parameter that would be consistent with the implied very high intrinsic luminosity of the source in this state. Thus, it is unlikely that the soft component needed to model epochs 1 and 2 in V404 Cyg comes from an accretion disk. In addition, V404 Cyg epoch 1 corresponds to the hardest X-ray state in Cyg X–3 (e.g., Koljonen et al. 2010). Cyg X–3 is a persistent wind-accreting source and thus always accreting at (fairly) similar mass accretion rates. The source likely stays in an intermediate state all the time, with an average power-law spectral index Γ ≥ 2.0. Thus, the spectral variation is at least partly due to changes in the surrounding medium. It has been suggested that the jet pressure can play a role in reducing the wind density in the line of sight, producing drastic variability in the spectral evolution of Cyg X–3 (Koljonen et al. 2018). However, in Cyg X–3, there is further evidence of a very soft state (hypersoft state), where the X-ray spectrum is dominated by a partially absorbed blackbody component with a very weak and flat hard X-ray tail (Koljonen et al. 2018). This might be (partly) similar to what is observed in V404 Cyg epoch 2, with the increased soft X-ray emission. The soft state in Cyg X–3 always precedes jet ejection episodes (Koljonen et al. 2010), which indicates that this might be the case for V404 Cyg as well. Koljonen et al. (2018) argued that in the hypersoft state, the jet turns off, allowing the stellar wind from the Wolf-Rayet companion to fill out the cavity created by the jet. This increases the density of the wind close to the X-ray source, providing a medium where multiple scattering can take place. When the jet is turned on later, it encounters a dense medium where an efficient energy dissipation can take place to boost the jet emission to Jy levels.

Gandhi et al. (2017) showed the multiwavelength evolution of V404 Cyg during the NuSTAR observation (their Fig. 1 and supplementary Fig. 4; their epoch 1 corresponds to our GTIs 10−14). They showed that during the change from the preflare state to the flaring state, there is a brightening in the 15 GHz flux at least by a factor of five and a change in the radio spectral index from negative values to above zero, indicating the start of the jet ejection. In model A, during epoch 2, there is a rise in the parameter evolution of R_in, Γ, log ξ, z, and inclination (Fig. 7). This indicates that the geometry of the reprocessor changes and moves farther away (the increase in inclination and R_in) with higher speeds (the increase in z), most likely caused by increased radiation pressure (the increase in log ξ) that also cools the electrons in the corona (the increase in Γ, and low values of kT_e). In model B, a similar evolution can be seen in z and in the blackbody flux (Fig. 8). Because the black body flux is proportional to the emitting area squared and temperature to the fourth (i.e., decreasing from epoch 1 to epoch 2; Table 7), this also implies an increase in the size of the blackbody emitter in addition to the increase in the speed of the scattering medium (even more than for model A). This means that what we may be seeing here is a jet ejection event following an accretion event that starts by pushing the reprocessing medium farther out and at the same changing its geometry from sphere-like to disk-like (Fig. 3) and exposing the inner accretion disk seen in the following high-luminosity flaring period (Walton et al. 2017; Gandhi et al. 2017). A similar spectral sequence can be seen as leading to a high-luminosity flare as well (see Walton et al. 2017; their Fig. 14), beginning from the Γ ∼ 2 low-cutoff spectrum with high absorption and evolving to a harder spectrum with a higher cutoff and low absorption (see Walton et al. 2017; their Fig. 15), but with a much faster evolutionary time (tens of seconds compared to ∼10 000 s, as shown in Fig. 2).

The intrinsic luminosity in epochs 1–2 is very high in model B, with an Eddington fraction of 0.2–0.3 in the 3–79 keV band, and when the model is extrapolated to lower energies, the model luminosity reaches the Eddington limit roughly in the 0.2–100 keV band. In the case of model A, we see only the scattered emission, while the intrinsic emission would be completely absorbed or scattered, and thus there is no reliable way to estimate the intrinsic emission. Therefore we assume that model B gives a better indication of the intrinsic luminosity. Based on the evolution of the model parameters outlined above, the super-Eddington luminosities, and the radio evolution, the likeliest scenario for V404 Cyg epochs 1–2 therefore is a super-Eddington accretion rate event that resulted in a large scale-height accretion flow and a powerful optically thick accretion disk wind that launched with mildly relativistic speed and led to a jet ejection event.

4.2. Is GRS 1915+105 in the hard state?

In July 2018, GRS 1915+105 entered an extended unusually low-flux X-ray phase followed by a change to a state with even lower average X-ray fluxes that were not seen before during the 27-year-long outburst (Fig. 13; light blue squares) but presented renewed flaring activity in radio as well as X-rays (Fig. 13; red diamonds and dashed lines). After the change to this peculiar state, radio monitoring data showed significant radio flaring (Motta et al. 2019) that started approximately at the time of the renewed X-ray activity. The radio flaring has since continued and is still ongoing at the time of writing (Trushkin et al. 2020). While this radio behavior is consistent with what has frequently been observed in the past (e.g., Fender et al. 1999; Punsly & Rodriguez 2013), this is the first time that significant radio activity does not seem to be associated with a strong X-ray counterpart.

Fig. 13. 2–20 keV daily light curve of GRS 1915+105 from The Monitor of All-sky X-ray Image/Gas Slit Camera (MAXI/GSC) since February 2019, with NuSTAR observations marked as arrows and radio flare detections (Motta et al. 2019; Trushkin et al. 2019; Koljonen et al. 2019) as vertical dotted lines. The data are colored and marked according to the spectral hardness shown in Fig. 14.

It is therefore not entirely clear whether the outburst of GRS 1915+105 is nearing its end or if the source is just highly obscured. The prolonged low-luminosity state since July 2018 does indicate that continuous super-Eddington accretion and Compton-thick wind production is unlikely. There seems to be no X-ray flare during the time of the flux drop (only later), which means that super-Eddington accretion and subsequent mass expulsion did probably not take place and did not cover the source. In addition, the X-ray spectrum is harder than ever observed from the source (blue points in Fig. 14), indicating that the source might have reached a regular hard state on its way toward quiescence.

Fig. 14. MAXI/GSC hardness-intensity diagram of GRS 1915+105 from daily monitoring observations since August 2009. The blue data points (dark triangles and light blue squares) indicate the recent low-luminosity state with increased spectral hardness. The light blue squares correspond to the anomalous state with occasional strong X-ray flares (red diamonds) and highly variable radio emission. The numbered green boxes correspond to the NuSTAR epochs.

On the other hand, the radio emission from the source has become more variable and presents flux densities that are among the strongest flares ever observed from GRS 1915+105 (Trushkin et al. 2020). In addition, sporadic X-ray flares display softer spectra and reach similar X-ray hardnesses as before the anomalous low-luminosity state (red points in Figs. 13 and 14), which at least in principle is consistent with varying absorption. While more detailed studies should be made in terms of radio to X-ray correlation, it seems that the X-ray emission remains at a rather constant level (Fig. 13), while large amplitude variations are evident in the radio monitoring data (Trushkin et al. 2019, 2020). The enhanced radio luminosity and variability can arise from additional energy dissipation either in merging shocks in the jet (as has been suggested as an explanation of the radio behavior of GRS 1915+105 in Vadawale et al. 2003) or the jet interacting with expelled matter from the accretion flow.

The unabsorbed luminosities of model A fits for V404 Cyg and GRS 1915+105 in the pre-flaring state are similar: ∼2 × 10³⁷ erg s⁻¹, indicating that a similar fraction of reflected emission is received from both sources. We have speculated that epochs 1–2 of V404 Cyg might be due to super-Eddington accretion, which might be then the case for GRS 1915+105 as well if the intrinsic emission is completely obscured to us. However, model B fits, where the absorbed line-of-sight component is taken into account, give a similar unabsorbed luminosity for GRS 1915+105 as model A fits, while for V404 Cyg, it is an order of magnitude higher. The resulting Eddington luminosity of 1–2% for both models of GRS 1915+105 is consistent with a regular XRB hard-state luminosity. In addition, the cospectrum for GRS 1915+105 in epoch 1 is markedly different from the latter epochs (or the cospectra of V404 Cyg and Cyg X–3; Fig. 12), displaying a low-frequency QPO and spectral power at least up to 10 Hz, while in others, there appears to be no power above 0.1 Hz, indicating that we have at least a partial view to the accretion flow in epoch 1. This is also consistent with the amount of the leaked or intrinsic emission received in the model B fits. It therefore seems more likely that GRS 1915+105 has reached a genuine hard X-ray state. However, it is clear from the X-ray spectra and the model fits above that the source is (partly) obscured in all NuSTAR epochs and likely continuously after MJD 58610. This may indicate that the accretion flow has changed from geometrically thin to thick, and that due to high inclination, it blocks the view to the central parts of the flow.

4.3. Implications for other sources

We have argued that the intrinsic X-ray emission from V404 Cyg, Cyg X–3, V4641 Sgr, and GRS 1915+105 is significantly affected by the surrounding medium in scattering processes. While it is clear that these four sources are unique among the XRB population with large accretion disks and high inclination angles, there is evidence that similar scattering takes place in other sources as well: (i) Strong and variable absorption has been found in the X-ray spectrum of Swift J1858.6−0814 during X-ray flaring that bears similarity to the flaring spectra from V404 Cyg and V4641 Sgr (Hare et al. 2020). In addition, P Cygni profiles were observed in the optical spectra, indicating a high-velocity wind (Munoz-Darias et al. 2019). (ii) The well-known accretion disk wind source GRO 1655–40 displays similar X-ray spectra as Cyg X–3 (Uttley & Klein-Wolt 2015). (iii) The hard X-ray emission from SS 433 has been speculated to be heavily reprocessed intrinsic X-ray emission from a supercritical accretion flow viewed through the optically thick accretion wind cone (Middleton et al. 2018). (iv) Other super-Eddington accretors, such as extragalactic ULXs, typically have soft X-ray spectra that might be similar to the epoch 2 spectrum of V404 Cyg. (v) Other high-mass XRBs, such as Cyg X–1, might exhibit reprocessing in the companion wind, and the dense environment around XRB jets can affect the properties of the jet emission through shocking, thus increasing the radiative efficiency and enhancing the radio luminosity.

Recently, puzzling observations of a low-luminosity soft state (< 0.01 L/L_Edd) in XRBs have been reported, including V4641 Sgr (Pahari et al. 2015), Swift J1753.5–0127 (Shaw et al. 2016), and 4UX-11630–47 (Tomsick et al. 2014). Pahari et al. (2015) showed that during a low-luminosity outburst or a renewed X-ray activity period in January–February 2014 (Tachibana et al. 2014; Uemura et al. 2014), the X-ray spectrum presented a soft state with reflection features including the ionized iron line and iron edges similar to the spectra presented in this paper. Swift J1753.5–0127 and 4U 1630–47 exhibited a low-luminosity soft state at the end of their outbursts in March–May 2015 and July 2010, respectively. With the low luminosity, it is difficult to attribute this component to a regular soft-state disk. The system parameters for these sources are not all well constrained, but there is some evidence for a high orbital inclination from optical variability of Swift J1753.5–0127 (Neustroev et al. 2014) and dipping phenomena in 4U 1630–47 (Kuulkers et al. 1998; Tomsick et al. 1998), while the orbital inclination of V4641 Sgr is fairly well constrained at 72° ± 4° (MacDonald et al. 2014). When we assume a high inclination, a significant change in the scale height of the accretion flow, for example, by disk warping or precession, might intercept our line of sight, resulting in strong absorption and reprocessing of the intrinsic emission that could thermalize the intrinsically hard spectrum. A misalignment between the black hole and accretion disk spins has been suspected for V4641 Sgr (Maccarone 2002; Gallo et al. 2014). This might cause the disk precession and perhaps explain the semi-regular interval of the weak outbursts (Negoro et al. 2018).

5. Conclusions

We have studied the NuSTAR and RXTE/PCA spectra of four unique XRBs, V404 Cyg, Cyg X–3, V4641 Sgr, and GRS 1915+105, which are known to present complex spectral evolution distinct from other XRBs. We showed that all sources have similar X-ray spectra at certain times that can be modeled by assuming that a Compton-thick medium surrounds the central X-ray source. This assumption is further enhanced by the fact that Cyg X–3 orbits its Wolf-Rayet companion star inside a high-density stellar wind that strongly affects the X-ray spectra. While the companion stars of V404 Cyg, V4641 Sgr, and GRS 1915+105 are not expected to present a high-density stellar wind, we attribute the obscuring medium to either a large scale-height accretion flow or to an optically thick equatorial outflow or envelope.

The results from fitting two physically motivated scattering models suggest that a low-luminosity phase preceding a flaring episode in the 2015 outburst of V404 Cyg is a heavily obscured, but intrinsically very bright (super-Eddington) accretion state. In this state, a dense medium fully covers the X-ray source, and the majority of the received flux comes from heavily absorbed intrinsic emission. After the fully obscured phase, the geometry changes to resemble a disk-like wind, and the majority of the received emission comes from the scattered or reflected component. During this time, large-amplitude flares are observed in X-rays and radio, indicating that some of the obscuring material has been removed either by a change in the accretion flow geometry, in the accretion wind geometry, or by the jet pressure. A part of the elevated emission may be due to the interaction of the jet with the obscuring matter. The shift of the iron line energy below 6.4 keV suggests that the scattering medium is in motion.

Similar spectral evolution to that of V404 Cyg is observed from the unusual low-luminosity state of GRS 1915+105, with the difference that the unabsorbed luminosity remains at a few percent of the Eddington luminosity. It is therefore more likely that the source has declined in flux and reached a regular hard X-ray state. Along with the state transition, the accretion flow must have thickened because there is evidence that the source is absorbed. Similarly, the weaker 2002 and 2003 outbursts of V4641 Sgr present similar spectra although the unabsorbed luminosities are very low, lower than 1% of the Eddington luminosity. In these cases, we might be seeing the sources through a disk wind or a geometrically thick accretion flow. Thus, this work highlights the importance of taking the reprocessing of the X-ray emission in the surrounding medium in the modeling of the X-ray spectra into account, which may well take place in multiple sources.

Acknowledgments

We thank Sara Motta for enlightening discussions and the anonymous referee for useful comments. KIIK was supported by the Academy of Finland project 320085. JAT acknowledges partial support from NASA under grant 80NSSC18K0574. This research has made use of data and software provided by HEASARC, which is a service of the Astrophysics Science Division at NASA/GSFC. The MAXI/GSC data has been provided by RIKEN, JAXA and the MAXI team.

References

All Tables

All Figures

Fig. 1. Collection of X-ray spectra from the four sources ordered into two panels according to their spectral shape. Top: NuSTAR spectra (FPMA/FPMB data both included) of V404 Cyg during the outburst period of 2015 (preceding intense X-ray flaring period), Cyg X–3 in the hard state, GRS 1915+105 during the anomalous low-luminosity state of 2019–2020 (preceding a period of sporadic X-ray and radio flaring), and the RXTE/PCA spectrum of V4641 Sgr during the outburst period of 2002. Bottom: NuSTAR spectra of V404 Cyg during the outburst period of 2015 (in between intense X-ray flaring), RXTE/PCA spectrum of V4641 Sgr during the outburst period of 2003, and NuSTAR spectra of GRS 1915+105 from the anomalous low-luminosity state (after the intense X-ray and radio flaring period). Note the spectral similarity between sources. The spectra of GRS 1915+105 and V4641 Sgr are renormalized by the amount shown for illustrative purposes.

In the text

Fig. 2. Initial modeling of the V404 Cyg preflare spectrum. Left: best-fit models fit to the hard X-ray data (10–79 keV), but plotted in the full data range. Different models are labeled, and the parameters are tabulated in Table 3. Solid lines refer to pure continuum models, dashed lines to absorbed continuum models, and dotted lines to reprocessed continuum models. Middle: fitting an absorbed cutoff power-law model to the full data range. The dot-dashed blue line corresponds to model A2 fit to the hard X-ray data (see left panel), the dotted green line corresponds to absorption and disk blackbody components added to the model, the dashed yellow line corresponds to partial absorption and smeared edge components added to the model, and the solid red line shows emission and absorption lines added to the model (parameters of the final model are tabulated in Table 4). See the text for more details. Right: fitting a reprocessed thermal Comptonization model to the full data range. The dot-dashed blue line corresponds to model R3 fit to the hard X-ray data (see the left panel), the dotted green line shows an increased value of the ionization parameter, the dashed yellow line corresponds to absorption and smeared edge components added to the model, and the solid red line shows partial absorption and an absorption line added to the model (parameters of the final model are tabulated in Table 4). See the text for more details.

In the text

Fig. 3. Geometry of the torus-shaped obscuring matter around the intrinsic X-ray source according to the parameter θtor in model B fits for all sources, in addition to a graphical representation of the column densities (gray, darker for a higher column density), the wind speeds (arrows), and the lines of sight (dashed lines). The thick dashed line around V404 Cyg epoch 1 and 2 corresponds to the partially absorbed blackbody component. The model parameters for each source can be found in Tables 6 and 8.

In the text

Fig. 4. NuSTAR 3–79 keV light curve of V404 Cyg divided into GTIs that total 30 000 counts (numbered). The numbers for GTIs that include count rates exceeding 1050 cts s−1 are not shown (excluding GTIs 13 and 14), and the data points are shown in gray with a reduced point size. The coloring scheme corresponds to the spectral shape shown in Fig. 5, with every other GTI shown in a different hue for clarity.

In the text

	Fig. 5. Changes in the spectral shape of V404 Cyg during the GTIs shown in Fig. 4. All the spectra have been normalized to the 7 keV flux in the first spectrum to show the difference in the spectral shape.
In the text

	Fig. 6. Spectra of V404 Cyg from GTIs 6–9 in Fig. 4, showing the change from a harder (6) to a softer spectrum (7–9) preceding the flaring period. The spectra are not normalized, thus showing the variable regimes below 6 keV and above 10 keV.
In the text

	Fig. 7. Model A parameter values for V404 Cyg spectra (individual GTIs in black, epochs in color).
In the text

	Fig. 8. Model B parameter values for V404 Cyg spectra (individual GTIs in black, epochs in color).
In the text

Fig. 9. Averaged NuSTAR FPMA data from V404 Cyg in the preflare (epochs 1–2) and outburst (epochs 3–5) stages as shown in Fig. 5, together with absorbed (top row) and intrinsic (bottom row) model B components. The total model (solid red line) consists of a sum of a blackbody component (dashed green line), a cutoff power-law component absorbed in and scattered off from the material in the line of sight (dotted blue line), scattered into the line of sight by the surrounding medium (dot-dashed light blue line), and unabsorbed or direct component (long-dashed yellow line). The middle panels show the residuals of the models to the data.

In the text

Fig. 10. Changes in the RXTE/PCA spectral shape of V4641 Sgr during the outburst periods of 1999, 2002, 2003, and 2005. All the spectra have been normalized to the 7 keV flux in the first spectrum to show the difference in spectral shape. Top panel: outburst spectra reminiscent of the spectra with high count rate and low absorption of V404 Cyg (Walton et al. 2017), bottom panel: spectra reminiscent of the flaring spectra with lower count rate and higher absorption of V404 Cyg in Fig. 5.

In the text

Fig. 11. Averaged NuSTAR FPMA data from GRS 1915+105 from the three epochs together with absorbed (top row) and intrinsic (bottom row) model B components. The total model (solid red line) consists of the sum of a cutoff power-law component absorbed in and scattered off from the material in the line of sight (dotted blue line), scattered into the line of sight by the surrounding medium (dot-dashed light blue line), and unabsorbed or direct component (long-dashed yellow line). Middle panels: residuals of the models to the data.

In the text

	Fig. 12. Cospectra of V404 Cyg, Cyg X–3, and GRS 1915+105 for NuSTAR epochs. The soft- and hard-band cospectra are shown separately when they display any differences. The best-fit single or broken power-law model is also plotted for all data. Additional Lorenzian functions are needed for the QPO and harmonics in the GRS 1915+105 epoch 1 data.
In the text

	Fig. 13. 2–20 keV daily light curve of GRS 1915+105 from The Monitor of All-sky X-ray Image/Gas Slit Camera (MAXI/GSC) since February 2019, with NuSTAR observations marked as arrows and radio flare detections (Motta et al. 2019; Trushkin et al. 2019; Koljonen et al. 2019) as vertical dotted lines. The data are colored and marked according to the spectral hardness shown in Fig. 14.
In the text

Fig. 14. MAXI/GSC hardness-intensity diagram of GRS 1915+105 from daily monitoring observations since August 2009. The blue data points (dark triangles and light blue squares) indicate the recent low-luminosity state with increased spectral hardness. The light blue squares correspond to the anomalous state with occasional strong X-ray flares (red diamonds) and highly variable radio emission. The numbered green boxes correspond to the NuSTAR epochs.

In the text