1 Introduction

The direct images of supermassive black holes in M87* and Sgr A* released by the Event Horizon Telescope revealed ring structures which are interpreted as advection dominated relativistically hot accretion flows in an immediate vicinity of the black hole event horizons (Event Horizon Telescope Collaboration et al. 2019, 2022). These images allow us to test complex accretion flow models with unprecedented precision. It is reasonable to ask whether a black hole accretion processes are universal and whether models constructed mostly for the supermassive black holes apply to accreting stellar-mass holes found in our Galaxy.

Stellar-mass black holes accreting matter from a stellar companion (hereafter XRB) are known to display several distinct X-ray spectral states characterized by different luminosities, X-ray spectral indices, multiwavelength components and timing properties (Fender and Belloni 2012; Belloni and Motta 2016). The two main XRB spectral states are a low-hard and a high-soft states. In the low-hard state, the luminosity is lower, the X-ray spectrum has a power-law shape and the source has a persistent radio component. During the high-soft state, the luminosity rises closer the Eddington limit, the spectrum is dominated by a thermal component and the radio counterpart is quenched. Sources can be also found in an intermediate state with thermal and power-law components both visible in the X-ray spectrum. Since many years it is generally understood that the X-ray soft component is produced by a thermal disk and the power-law component is emitted by the so called “corona” which has some, not yet well understood, connection to a relativistic jet emitting photons in the radio (Done et al. 2007). One of the long-standing questions is what is the geometry and properties of this corona and whether the corona is actually a compact jet or rather a part of the disk and, if so, then how is the corona connected to the jet. All these questions are difficult to answer solely observationally because accretion flows in XRBs are spatially unresolved.

In this work we construct a more detailed model of the accretion flow in XRB in the low-hard state using state-of-the-art general relativistic radiation-magnetohydrodynamic (GRRMHD) simulations. Recent Imaging X-ray Polarimetry Explorer (IXPE) measurements of polarization of X-ray emission in XRBs in a hard state (e.g., Krawczynski et al. 2022) opened up a unique opportunity to test accretion flow scenarios in these systems as suggested already by Lightman and Shapiro (1976) and Rees (1975). In this work we show that the GRRMHD simulations indicate consistency with these new observational constraints.

The paper is structured as follows. In Sect. 2 we describe our model for hard-state accretion. In Sect. 3 we show the polarimetric properties of the model. We compare our modeling results to observations and to other models in Sect. 4. We discuss prospects and directions for future studies in Sect. 5.

2 Accretion model for a hard spectral state

We carry out GRRMHD simulations of low-hard state accretion flow onto stellar mass black hole using first-principles GRRMHD code ebhlight (Ryan et al. 2017). The state-of-the-art code couples evolution of plasma and magnetic fields to the evolution of radiation field accounting for radiative looses and radiative forces. The radiation model includes continuum processes such as: synchrotron emission, bremsstrahlung emission, self-absorptions and Compton scatterings of photons. Additionally, the code follows ions and electron temperatures separately allowing for development of two-temperature plasma flows (Ressler et al. 2015). The simulations include evolution of radiation field and hence they are not scale free. Inspired by recently published results on X-ray polarization in XRB Cygnus X-1 (Krawczynski et al. 2022), we scale our model exactly to this system. In our simulation we set black hole mass to \(M_{\mathrm{BH}}=20\mathrm{M_{\odot }}\) (Miller-Jones et al. 2021) and dimensionless spin of the black hole to \(a_{*}=0.9375\) (Fabian et al. 2012; Gou et al. 2014). The GRRMHD simulations require specifying the gas density and magnetic field strength via scaling parameter \(\mathcal{M}_{\mathrm{unit}}\). We run two simulations with \(\mathcal{M}_{\mathrm{unit}} = 1.5 \times 10^{11}\) and \(2.5 \times 10^{11}\). Two values of \(\mathcal{M}_{\mathrm{unit}}\) are chosen to model emission for observers with low (disk face-on) and high (disk edge-on) viewing angles. The mass accretion rates in both models are at the level of \(\dot{m}\equiv \dot{M}/\dot{M}_{\mathrm{Edd}}\sim {\mathrm{few}} \times 10^{-4}\). In X-rays both models are as bright as \(L_{X} \approx 2 \times 10^{36}\) ergs/s consistent with luminosity of Cygnus X-1 during IXPE detection reported by Krawczynski et al. (2022).Footnote 1 The Bolometric luminosities of our models are \(L_{\mathrm{Bol}} \approx 10^{38} \,{\mathrm{ergs\,s^{-1}}}\) (\(L_{\mathrm{Bol}}/L_{ \mathrm{Edd}}\approx 0.05\)).

The simulation initial conditions assume torus in hydrodynamical equilibrium in the equatorial plane of the black hole. The torus is seeded with a loop of poloidal magnetic field. Further details of the simulation setup are provided in the Appendix A.1. It is important to note here that our models initially do not contain a truncated geometrically thin, optically thick diskFootnote 2 component (Shakura and Sunyaev 1973). This means that in our models there is no reflection of the hot accretion flow from the thermal thin disk and hence there is no thermal photons scattering off the hot corona. Our simulations self-consistently evolve magnetic fields. The magnetorotational instability leads to the development of a turbulent accretion flow. The accretion flow is accompanied by wind and a jet with organized magnetic fields. The accretion flow models studied here are magnetically arrested (MADs, Tchekhovskoy et al. 2011, McKinney et al. 2012, and see Appendix A.1). The powers of a jets produced by the MAD simulations are consistent with the estimated power of the relativistic jet observed in Cyg X-1 (\(P_{ \mathrm{jet}}=10^{36}-10^{37}\) ergs/s, Gallo et al. 2005, Russell et al. 2007).

Figure 1 presents maps of plasma density, electron temperature, plasma \(\beta \) (gas to magnetic pressure) and radiation \(\beta _{r}\) (gas to radiation pressure) of the very inner accretion flow at a later time moment of evolution. The very inner accretion flow is geometrically thick, it has two-temperature structure and it is rather strongly magnetised. The optical thickness of the simulations is at most \(\tau _{es} \approx n_{e} \kappa _{es} r_{g} \approx 0.03\) hence the radiation emitted from the gas may have impact on the gas via radiative back reaction rather than Compton forces. Interestingly, in the disk body the regions of colder electrons are spatially coincidental with low plasma \(\beta \) regions suggesting that electrons are cooling rapidly in regions of stronger magnetic fields. These sub-relativistic electrons may enhance Faraday rotation of polarized synchrotron photons (seed photons for Compton process) which may impact the polarization properties of the X-ray emission.

Fig. 1
figure 1

Geometry of the inner accretion flow in the hard spectral state in GRRMHD simulation with \({\mathcal {M}}_{\mathrm{unit}}=1.5 \times 10^{11}\) at time moment \(t=4000\) M. Panels from left to right show the plasma density (in code units), the dimensionless electron temperature, the plasma \(\beta \) parameter indicating dynamically important magnetic fields, and the radiation \(\beta _{r}\) parameter (ratio of gas pressure to radiation pressure). The black contours follow the magnetic field lines

3 Mock spectral energy distribution and polarization

Although the GRRMHD simulations do produce self-consistent emission spectra those spectra do not have polarization information. Hence, we post-process our GRRMHD models with radpol, recently developed, fully relativistic, Monte Carlo code for polarized radiative transfer (Mościbrodzka 2020, 2022). The details of the calculation are provided in the Appendix A.2.

Figure 2 shows polarized broadband SED of the GRRMHD models for a couple of viewing angles around times \(t=4000~{\mathrm{M}}\). The lower viewing angle model has slightly higher \({\mathcal {M}}_{\mathrm{unit}}=2.5\times 10^{11}\) to match the observed X-ray flux and spectral slope. At the shown time the simulation is already relaxed and approximately steady-state with variations due to magnetic turbulence. Each SED has three components: synchrotron hump, X-ray scattered emission and bremsstrahlung of which synchrotron hump and Compton humps are the most prominent. The bottom panels in Fig. 2 show polarimetric information. The fractional linear polarization \(|m|=\sqrt{Q^{2}+U^{2}}/I\) in X-ray band is small, less than 10% for both viewing angles. The fractional linear polarization of the synchrotron emission is naturally higher. The electric vector position angle (\(EVPA \equiv 0.5arg(Q+iU)\)) is zero for close to edge-on inclination angles, meaning that EVPA is aligned with the jet for nearly edge-on view of the model. For lower viewing angles, the EVPA significantly deviates from zero and displays some frequency dependency.

Fig. 2
figure 2

SEDs of the GRRMHD models at a viewing angle for which the SED matches observational constraints (\(i=90\) deg, left panel) and at a viewing angle equal to the estimated inclination angle of the binary system orbital plane (\(i=30\) deg, right panel). Top panels show total intensity spectra, middle panels show the amplitude of the fractional linear polarization \(|m|\%\) and lower panels display the linear polarization angle, EVPA. The gray shaded region marks bandwidth of the X-ray observations (2-8 keV)

Unfortunately, it is not straightforward to measure how much of the X-ray polarization observed in the model is caused by the scatterings or the seed photon polarization because one cannot disentangle different effects easily in these complex simulations (e.g., the scattering process does depend on the polarization of incident photons). Due to the complexity of the model it is even more difficult to explain the behaviour of the EVPA as a function of viewing angle and frequency.

In addition to the spectrum our radiative transfer simulations track the origin of emission detected in different spectral windows. In Fig. 3 we show a snapshot of the GRRMHD simulation together with emission maps produced by our radpol scheme. Notice that in comparison to Fig. 1, Fig. 3 shows the model on larger scales (within 50 M) and plasma density and magnetic field strength are shown in c.g.s. units. Here it is evident that the electrons in the GRRMHD model are mildly relativistic with temperature increasing with latitude. Hence the seed synchrotron emission is produced in the inner parts of the accretion disk but mostly along the jet funnel wall. The inverse-Compton scatterings take place along the jet wall but also in the disk region.

Fig. 3
figure 3

Example maps of plasma density, magnetic field strength (both scaled to Cygnus X-1 system and given in c.g.s. units), dimensionless electron temperature and maps of origin of photons at different energy bands in one snapshot of the GRRMHD simulation. The maps of the photon origin are \(\phi \) angle averaged to improve signal-to-noise ratio in the maps. The color scale in the emission maps (two right-most panels) is proportional to the number of photons emitted. Notice that the size of the maps shown here is larger in comparison to Fig. 1

4 Comparison with observations and other models

The IXPE detected X-ray fractional polarization of Cygnus X-1 in hard state to be on average \(|m|=\)4% (slightly rising with frequency) and the EVPA to be aligned with the radio jet position angle (in our reference system this translates to EVPA=0). Given the simplifications in the model (e.g., the simulations are axisymmetric) our GRRMHD simulations of the hard state are quite consistent with these observational constraints. Low polarization fraction and EVPA aligned with the black hole spin or jet axis can be simultaneously recovered by GRRMHD simulations assuming nearly edge-on viewing angles. In our model the X-ray emission is emitted predominantly by the jet wall although some disk emission is also contributing. For viewing angle of \(i=30\) deg (estimated viewing angle of the Cyg X-1 orbital plane, Miller-Jones et al. 2021) the EVPA can become nearly perpendicular (EVPA\(\approx -100\) deg) compared to the observed value. Overall our results indicate a possible misalignment between the jet axis (black hole spin axis) and the angular momentum of the binary or the outer accretion disk.

In our model the X-ray emission and its polarization are produced by a partially outflowing matter (jet wall). This is consistent with ideas recently put forward by others (Poutanen et al. 2023; Dexter and Begelman 2023) however the first-principles simulations presented here prefer nearly edge-on viewing angles. Our results are in contradiction with the IXPE data interpretation presented in Krawczynski et al. (2022) where authors argue that the measurement made in Cygnus X-1 can be interpreted as a “sandwich”-like corona above the thin disk or “corona” alone (with either external Comptonization or self-synchrotron Comptonization) observed at intermediate inclinations angles (\(i=45-60\) deg). For such high viewing angles our model are depolarized in X-rays and EVPA deviates from zero. All these discrepancies can be attributed to different assumptions made when using different modeling techniques. For example, here we do not include the thin disk component, in contrast to the previous work. On the other hand, we include the previously neglected Faraday rotation which very likely plays a role in depolarization process of synchrotron emission in the disk where sub-relativistic electrons exist. Our models also include self-consistently evolved magnetic field geometry which is not included in the aforementioned geometrical models. Alternatively, the current model may be incorrect and another GRMHD scenario should be considered to explain the IXPE measurements.

5 Future prospects

Overall, the accurate detection of the X-ray polarimetry is a powerful new tool for constraining GRMHD simulations of XRBs and accreting supermassive black holes. In this work, we have made a first attempt to model polarimetric characteristics of XRB, utilizing two state-of-the-art numerical techniques. GRRMHD simulations, although already complex, still have several uncertainties that may impact the interpretation of the observations. At least two types of developments could be pursued.

The sub-grid electron thermodynamics in the GRMHD simulations is uncertain. One could consider alternative to currently used electron heating model which assumes purely Alfvénic cascade for energy dissipation. However it is not clear which model for electron heating is realistic for a given GRRMHD simulation (e.g., Howes 2010; Kawazura et al. 2018; Satapathy et al. 2023). Another missing component in the GRRMHD models are the non-thermal electrons which are likely accelerated along the jet. Including the effect of acceleration using first-principles would require developing new modeling approaches.

A more feasible goal would be to calculate broadband polarization of a fully three dimensional simulations such as the ones presented in Dexter et al. (2021). Our current model in an axisymmetric version of fully three dimensional models M3/M4 in Dexter et al. (2021) but notice that the our model differs in the black hole mass, evolution time and the size of the gas-radiation interaction zone (assumed to be 40 M by Dexter et al. (2021)). Three dimensional simulations allow to study emission from models with accretion disk tilted with respect to the equatorial plane so they may be helpful to explain the discrepancy between the viewing angle of the inner accretion flow and the orbital plane of the system. Finally, only three dimensional GRRMHD models are suitable for time variability studies.