New Constraints on the Spin of the Black Hole Cygnus X-1 and the Physical Properties of its Accretion Disk Corona

We present a new analysis of NuSTAR and Suzaku observations of the black hole Cygnus X-1 in the intermediate state. The analysis uses kerrC, a new model for analyzing spectral and spectropolarimetric X-ray observations of black holes. kerrC builds on a large library of simulated black holes in X-ray binaries. The model accounts for the X-ray emission from a geometrically thin, optically thick accretion disk, the propagation of the X-rays through the curved black hole spacetime, the reflection off the accretion disk, and the Comptonization of photons in coronae of different 3-D shapes and physical properties before and after the reflection. We present the results from using kerrC for the analysis of archival NuSTAR and Suzaku observations taken on May 27-28, 2015. The best wedge-shaped corona gives a better fit than the cone-shaped corona. Although we included cone-shaped coronae in the funnel regions above and below the black hole to resemble to some degree the common assumption of a compact lamppost corona hovering above and/or below the black hole, the fit chooses a very large version of this corona that makes it possible to Comptonize a sufficiently large fraction of the accretion disk photons to explain the observed power law emission. The analysis indicates a dimensionless black hole spin parameter a between 0.86 and 0.92. The kerrC model provides new insights about the radial distribution of the energy flux of returning and coronal emission irradiating the accretion disk. kerrC furthermore predicts small polarization fractions around 1% in the 2-8 keV energy range of the recently launched Imaging X-ray Polarimetry Explorer.


INTRODUCTION
The 2020-2030 decade promises to lead to several breakthrough discoveries in the field of astrophysical studies of black holes and black hole accretion. The pointed NuSTAR (Harrison et al. 2013), NGO (Gehrels et al. 2004), NICER (Gendreau et al. 2016), Chandra (Havey et al. 2019), and XMM-Newton (Kirsch et al. 2017) X-ray observatories will continue to operate, and will be joined by the X-ray and γray missions IXPE (Weisskopf et al. 2021), XL-Calibur (Abarr et al. 2021), and COSI (Siegert et al. 2020) with polarimetric capabilities, and by the XRISM mis-black holes are thought to accrete through a geometrically thin, optically thick accretion disk surrounded by hot coronal plasma (e.g., Remillard & McClintock 2006;McClintock et al. 2014). The accretion disk emits a multi-temperature thermal component peaking at keVenergies which follows roughly the analytical models of Shakura & Sunyaev (1973) ;Novikov & Thorne (1973). The Comptonization of the emission in a hot corona is believed to explain the non-thermal continuum emission extending from a few keV to a few ten keV or even a few hundred keV. Returning and coronal emission irradiating the disk gives rise to a reflection component, including the kinematically and gravitationally broadened fluorescent Fe K-α emission around 6.4 keV (Fabian et al. 1989). Some X-ray spectroscopic observations show furthermore evidence for winds (Miller et al. 2016), and some hard X-ray polarization and γ-ray observations suggest the presence of a collimated plasma outflow (jet) (e.g., Zdziarski et al. 2017).
We focus in the following on the disk, the coronal, and the reflected emission. As these emission components originate close to the black hole, they can inform us about the extremely curved spacetime close to the event horizon of the black hole, and allow us to constrain the inclination (the angle between the spin axis and the observer) and the spin of the black hole (McClintock et al. 2014;Miller 2006). As the fidelity of our models continues to improve, test of general relativity based on Xray observations may become feasible (Johannsen 2016;Krawczynski 2018;Bambi et al. 2021). As the accretion flow dissipates most of its energy close to the black hole, the three emission components are furthermore excellent diagnostics of the accretion flow physics. They present us with the opportunity to test models of the vertical structure of the accretion disk, the angular momentum transport inside and outside of the disk, the transfer of energy between ions, electrons, magnetic field, and radiation, and the properties of the accreted plasma and magnetic field (e.g. Yuan & Narayan 2014).
Current state-of-the-art analyses (e.g. Tomsick et al. 2018) fit the energy spectra with a combination of a multi-temperature disk model such as KERRBB (Li et al. 2005) and DISKBB (Mitsuda et al. 1984;Makishima et al. 1986;Kubota et al. 1998) plus a powerlaw or Comptonized powerlaw (e.g., Zdziarski et al. 1996;Życki et al. 1999) and a model of the reflected emission. The energy spectra of the reflected emission are calculated with models such as reflionx, reflionx_hd (Ross & Fabian 2005;Ross et al. 1999), or XILLVER (García et al. 2011(García et al. , 2013(García et al. , 2014 that are based on radiative transfer calculations in a plane parallel atmosphere. The emission components are convolved with a kernel to account for the frequency shifts from the relativistic motion of the emitting and reflecting plasma and the gravitational frequency shift incurred by the Xrays as they climb out of the curved Kerr spacetime (e.g. Dauser et al. 2010).
The RELXILL model of Dauser et al. (2014); García et al. (2014) assumes a point source of hard X-rays positioned at a height h on the rotation axis of the black hole. The lamppost model predicts the dependence of the flux irradiating the accretion disk as a function of radial distance r, but cannot predict the absolute coronal flux, the energy spectrum of the coronal emission, nor the polarization of the coronal emission. The model assumes that the coronal energy spectrum and the density of the reflecting plasma follow simple parameterizations. The lamppost assumption makes it possible to account for the dependence of the energy spectrum of the reflected emission on the direction of the reflected photons.
kerrC replaces the lamppost hypothesis with spatially extended coronae of different shapes and locations, and different physical properties (electron temperature and density). Modeling the Comptonization of the accretion disk emission in the corona from first principles, kerrC can predict the absolute flux of the Comptonized emission, and the radial dependence of the flux and energy spectrum of the coronal emission irradiating the accretion disk. The modeling with kerrC can be used to measure the black hole spin and inclination and to constrain the shape and physical properties of the corona.
kerrC implements the following features: • It includes a wide range of pre-calculated models: we simulated black hole spin parameters a between -1 (maximally rotating black hole with counterrotating accretion disk) and +0.998 (near maximally spinning black hole with co-rotating accretion disk), for widely varying coronal geometries (wedge and cone-shaped coronae of different locations and sizes), and physical parameters (electron densities and temperatures).
• kerrC makes absolute rather than relative predictions of the thermal, Comptonized, and reflected emission. This allows for a comprehensive test of the underlying assumptions. For example, kerrC can be used to evaluate if a corona can be very close to the black hole as indicated by some spectral and timing studies (e.g. Wilms et al. 2001;Fabian et al. 2009;Chiang et al. 2015;Uttley et al. 2014), and, at the same time be sufficiently large to Comptonize enough accretion disk photons to account for the observed power law and reflected emission components (Dovčiak & Done 2016;Ursini et al. 2020). Furthermore, fitting the continuum and the lines at the same time, kerrC can constrain system parameters such as the black hole spin and inclination more tightly than models that use additional fitting parameters.
• The model tracks the thermal, coronal, and returning emission, and the emission reflected off the accretion disk.
• kerrc accounts for the Comptonization of the emission following the reflection off the disk (see the related discussion by Steiner et al. 2017).
• kerrC uses pre-calculated geodescis which are convolved with the chosen reflection model "on the fly" while fitting the data. The architecture makes it possible to convolve the kerrC configurations with the XILLVER models, and to introduce and fit additional parameters. kerrC allows for example to fit a parameter describing the dependence of the photospheric electron density on the distance from the black hole.
• kerrC calculates the flux (Stokes I) and polarization (Stokes Q and U ) energy spectra. The model can thus be used for the analysis and interpretation of the flux and polarization energy spectra acquired with the recently launched Imaging Xray Polarimetry Explorer (IXPE) mission and the upcoming X-ray polarization mission XL-Calibur.
The approach described here is complimentary to parallel efforts that combine first-principle GRMHD or 2tGRRMHD simulations with raytracing calculations (e.g., Kinch et al. 2016Kinch et al. , 2019Kinch et al. , 2020Kinch et al. , 2021Liska et al. 2018, 2019a,b, 2020. We envisage that the results from firstprinciple simulations can be used in future to refine the kerrC model. Conversely, the results from fitting observational data with kerrC can be used to identify shortcomings of the first-principle simulations. The interested reader is referred to (Sunyaev & Titarchuk 1985;Haardt et al. 1994;Nagirner & Poutanen 1994;Poutanen 1994;Poutanen & Svensson 1996;Zdziarski et al. 1996;Poutanen et al. 1997;Życki et al. 1999;Schnittman & Krolik 2010;Dauser et al. 2013;Gonzalez et al. 2017;Beheshtipour et al. 2017;Zhang et al. 2019) for earlier analytical and numerical studies of the properties of the coronal emission.
The rest of the paper is organized as follows. We describe the numerical simulations in Sect. 2, and the kerrC implementation in Sect 3. We present the results from fitting Suzaku and NuSTAR observations of Cyg X-1 in Sect. 4. We present predictions for IXPE in Sect. 5. We summarize and discuss the results in Sect. 6.
In the following we use c = 1, and define the gravitational radius of the black hole of mass M to be r g = GM .

Overall Architecture and Simulated Configurations
The kerrC fitting model is based on a library of 68,040 raytraced black hole, accretion flow, and corona configurations. For each configuration, 20 million events are generated. After describing the raytracing simulations in this section, we detail how the simulated data are used in the fitting code in Section 3.
The parameter grid describing the simulated configurations is summarized in Table 1. Black holes with spin parameters a between -1 and Thorne's theoretical maximum spin a = 0.998 (Thorne 1974) are simulated. The sampling is denser close to the maximum spin as the accretion disk properties change rapidly as a approaches 1. We perform the simulations for a black hole with a mass M of 10 M ⊙ , and with different temperature profiles. Taking advantage of the analytical result that the temperature profiles depend on the accretion rate only through a multiplicative scaling factor σ (Page & Thorne 1974), we simulate each model for ten σ values. The simulations with different temperature profiles can then be used to derive energy spectra for different black hole masses and accretion rates (see Sect. 3).
We simulate two families of coronal geometries. The first family are wedge-shaped coronae surrounding the accretion disk (Fig. 1). The coronae extend from the innermost stable circular orbit r 1 = r ISCO (Bardeen et al. 1972) to r 2 equal 25, 50 or 100 r g with a half-opening angle θ C between 5 • and 85 • above and below the accretion disk. Note that the wedge-shaped coronae morph into near spherical coronae as θ C increases from 5 • to 85 • .
The second family are cone-shaped coronae in the funnel regions around the black hole spin axis (Fig. 2). The cone-shaped coronae extend from r 1 = 2.5 r g , 10 r g , or 50 r g to r 2 20 r g , 50 r g , or 100 r g , respectively, above and bellow the black hole. The opening angles θ C range from 5 • to 45 • . For θ C = 5 • , the cone-shaped coronae resemble a jet in the funnel region above and below the black hole. For θ C = 45 • , the cone-shaped coronae resemble less collimated structures.
Corona electron temperatures T C between 5 keV and 500 keV, and integrated optical depths between the case of no corona (τ C = 0) and τ C = 2 are simulated. For wedge-shaped coronae, the optical depth is measured vertically from the accretion disk to the upper (or lower) edge of the corona. For the cone-shaped coronae, the optical depth is measured radially from the inner to the outer edge of the corona. For a fast spinning black hole (a = 0.998) with a wedge corona extending from the ISCO to 100 r g , the highest temperature (T C = 500 keV), an optical depth perpendicular to the accretion disk of (τ C = 2), and a wedge angle of 20 • , the code produces a power law component with a photon The code tracks photons forward in time, and allows us to collect photons at arbitrary locations of the observer relative to the accreting system. A single simulation set can thus be used to infer the observations for all black hole inclination angles.

Raytracing Code
We generate the simulated energy spectra with the code xTrack described in (Krawczynski 2012;Hoormann et al. 2016;Beheshtipour et al. 2017;Krawczynski et al. 2019;Abarr & Krawczynski 2020a. We give here a slightly updated description of the code that includes several recent improvements. The code uses the Kerr metric g µν in Boyer-Lindquist (BL) coordinates x µ = (t, r, θ, φ). Photons are emitted and scatter in the reference frames of the emitting or scattering plasma, and are transported in the global BL coordinate frame. The reference frame transforma-tions use orthonormal tetrads. The code assumes that the black hole accretes through a geometrically paper thin, optically thick accretion disk (Shakura & Sunyaev 1973;Novikov & Thorne 1973;Page & Thorne 1974). For each black hole spin we calculate a fiducial temperature profile T 0 (r) with the accretion rateṀ chosen such that the luminosity is 10% of the Eddington luminosity In the above equations, η is the fraction of the gravitational energy of the matter that is converted to heat as the matter accretes from infinity to r ISCO in units of the matter's rest mass energy. The fraction is given by: with E † (r ISCO ) being the specific energy at infinity of the matter at r = r ISCO : where g is the Kerr metric and u µ the four velocity of the matter at the ISCO. Given the angular velocity dφ/dt = Ω = ±M 1/2 /(r 3/2 ± aM 1/2 ) ( Bardeen et al. 1972), we get u µ = u t (1, 0, 0, Ω) from the condition u 2 = −1. The emissivity (power per comoving time and area) F (r) is then given by Equations (11)-(12) of (Page & Thorne 1974). We assume a fiducial temperature profile with  where σ SB is the Stefan Boltzmann constant. As mentioned above, we perform the simulations for a range of scaled temperature profiles: with σ between 0.05 and 2. We generate initial photon energies assuming a diluted blackbody en-ergy spectrum with a hardening factor of f h = 1.8 (Shimura & Takahara 1995).
For each accretion disk and corona configuration, we simulate 2×10 7 photons, using 10,000 logarithmically spaced radial bins with r running from r ISCO to 100 r g . The photons are launched into the upper hemisphere with constant probability per solid angle with a limb brightening weighting factor and initial polarization given by Table XXIV of (Chandrasekhar 1960). We normalize the limb brightening weighting factor so that the average weighting factor is 1.
An adaptive step-size Cash-Karp integrator based on a 5'th order Runge-Kutta algorithm is used to solve the geodesic equation, to update the the position (x µ ), and the photon's wavevector k µ , and to parallel transports the polarization vector f µ (Abarr & Krawczynski 2020a,b). The polarization fraction remains constant between scattering events. The step size is reduced if the traversed coronal optical depth evaluated in the rest frame of the corona exceeds 2% of the total optical depth of the corona. The step size is furthermore modified if the photon enters or leaves the corona, so that the end of the step coincides with the corona boundary.
Photons are tracked until their radial BL coordinate drops below 1.02 times the r-coordinate of the event horizon (at which point we assume that the photon will disappear into the black hole) or reach a fiducial observer at r obs = 10,000 r g . In the latter case, the wave and polarization vectors are transformed into the reference system of a coordinate stationary observer, and key information is written to disk.

Calculation of Absolute Fluxes Reaching the
Observer and Reaching the Accretion Disk Each event is counted with a weighing factor that allows us to make absolute flux predictions. As mentioned above, the thin disk solution of Page & Thorne (1974) gives the plasma-frame energy flux F (r) of the photons emitted at radius r. The number of photons emitted per Boyer-Lindquist dt, dr and dφ is given by (Krawczynski 2012): with −g trφ (r) is the square root of the negative of t−r −φ-part of the metric. For the Kerr metric, the factor simplifies to −g trφ (r) = r in the equatorial plane. The mean plasma frame energy of the thermally emitted photons is given by Although Equation (7) has been derived from a relativistic invariant, the interpretation is simple: the photon flux per unit coordinate time dN/dt equals the proper area-time volume element of the emitting disk segment −g trφ (r) dr dφ dt times the emitted energy flux per unit area F (r) divided by the product of the mean energy of the emitted photons < E(r) > times dt.
The i'th disk segment extending from r i − ∆r i /2 to r + ∆r i /2 thus emits photons at a rate of: (9) Given that we generate N bin events per radial bin, each simulated photon transports the photon rate dNi dt /N bin . The photon rate per simulated event launched from the i'th bin is thus give by the weight: For a source at the distance d from us, the flux spreads over an area: A Ω = ∆Ω d 2 (11) before reaching us. Here, ∆Ω is the solid angle of the θ-window that we use to collect the simulated photons. Each simulated event thus transports the photon flux: per unit observer time and area. Differential photon fluxes per unit time, unit area, and unit energy (Stokes I, Q, and U ) can be obtained by binning the results in observer energy and dividing the sum of all f i in each bin by the width of the energy bin. The calculation of the reflected energy spectrum requires the knowledge of the flux and the photon index of the emission irradiating the j'th disk element. Tracing the argument that led to Equation (10) backwards, we infer that each event being launched from the i'th bin of width ∆r i that reaches the j'th bin of width ∆r j adds a plasma frame energy flux of to the energy flux irradiating the bin. Here, E d is the plasma frame energy of the photon irradiating the disk. This equation has again a simple interpretation: the energy flux transported by each event starting at bin i and reaching bin j is the energy flux F (r i ) leaving bin i divided by the number of events launched from the i'th bin times the ratio of the area-time volumes of the i'th and j'th radial bins times the fractional energy gain or loss of the photon between emission from bin i and arrival at bin j.
Summing f i→j over all events irradiating the j'th radial bin gives the total plasma frame energy flux F x (r j ) irradiating the bin. Reflections off the disk and/or Compton scattering processes modify the statistical weight of a photon.
In a photospheric plasma with comoving electron density n e (Ross & Fabian 1993;Ross et al. 1996;Ross & Fabian 1996;Ross et al. 1999), the ionization can be characterized with the ionization parameter ξ. The parameter is proportional to the ratio of the photoionization rate (∝ n e ) divided by the recombination rate (∝ n 2 e ): The energy spectra of all photons arriving in the radial bin j can be used to fit the photon index Γ j (from dN/dE ∝ E −Γ ) of the emission irradiating the bin.
A limitation of our code should be mentioned: photons reflecting off the disk more than once are not modeled fully self-consistently. For these photons, only their first encounter with the disk contributes to the calculation of ξ and Γ, and additional encounters are neglected. The statistical weight of the subsequent photon-disk encounters depends on the ξ and Γ values used for the previous encounters. Properly accounting for multiple diskphoton interactions would thus render the problem nonlinear. For the purpose of calculating the predicted energy spectra, kerrC does track the photons through multiple disk encounters. However, the XILLVER tables are only used to modify their statistical weight for their first disk encounter. For subsequent encounters, the weights from the Chandrasekhar treatment are used. We estimated the impact of these approximations on the fitted parameters, by using the "opposite" assumptions, i.e. by including the additional photon-disk encounters in the ξ and Γ calculation with their full Chandrasekhar weight, and by completely excluding photons scattering more than once from entering the predicted energy spectra. We find that the approximations have a negligible impact in case of the cone-shaped corona models as few photons return multiple times to the disk. Adopting the alternative assumptions for the wedge-shaped coronae tends to harden the energy spectra impinging on the disk for the high-ξ models and to soften the predicted energy spectra for the low-ξ models. We will explore the problem and possible solutions in greater detail in future papers.

Comptonization
We use the Comptonization code of (Beheshtipour et al. 2017;Beheshtipour 2018). The coronal plasma is assumed to be stationary in the Zero Angular Momentum Observer (ZAMO) frame (Bardeen et al. 1972). See (Krawczynski 2021) for a formalism to implement relativistically moving coronal gas. For each integration step, the traversed Thomson optical depth ∆τ is calculated and a random number generator decides if the photon is considered for a coronal scattering event.
The scattering event is simulated by transforming the photon's wave and polarization four vectors into the ZAMO frame, and subsequently into the frame of a scattering electron. We assume that electrons move isotropically in the ZAMO frame with energies drawn from a relativistic Maxwell-Jüttner distribution (Jüttner 1911). The scattered photon direction is drawn from a distribution with equal scattering probability in the electron rest frame. Subsequently, the Stokes parameters are calculated relative to the scattering plane. The fully relativistic Fano scattering matrix derived from the Klein-Nishina cross section (Fano 1949;McMaster 1961) is used to calculate the Stokes vector of the outgoing beam. We use the Stokes I component of the scattered beam and a random number generator to decide if a photon actually scatters. This rejection technique enables us to account for the dependence of the scattering probability on the scattering direction and on the photon energy. If the photon indeed scatters, the statistical weight is multiplied with the kinematic factor 1 − cos (θ) with θ being the pitch angle of the incoming photon and electron, accounting for the relative motion of the photons and electrons (Beheshtipour et al. 2017). The scattering is followed by the back transformation of the photon wave and polarization vectors into the BL frame.

Preliminary Disk Reflection
For photons irradiating the accretion disk, we perform a simple reflection in the reference frame of the disk plasma based on the analytical results of Chandrasekhar for the reflection of a polarized beam of photons off an indefinitely deep plane-parallel atmosphere (Chandrasekhar 1960, Section 70.3, Equation (164) and Table XXV). The reflection changes the statistical weight of the reflected beam, and the linear polarization fraction and angle. For each reflected photon, the code saves information about the incident photon and the reflected photon as measured in the disk reference frame, and information about the photon when it reaches the observer. This information is used during the actual fitting of the observational data to re-weigh the reflected photons, so that the reflected energy spectrum resembles that from the XILLVER radiation transport calculations for the self-consistently derived ionization parameter ξ(r) and the spectral index Γ(r) of the emission irradiating the accretion disk.

Implementation Details
The kerrC fitting model uses a data bank of 1.36 trillion simulated events (20 million events for each of the 68,040 configurations.) The raytracing code stores for this purpose three types of data: data of photons reaching the observer without reflecting off the disk (direct emission data), data of photons irradiating the disk (disk data), and data of photons reaching the observer after scattering at least once off the accretion disk (reflected emission data). This classification is independent of the scattering of photons in the corona. The data are stored in Hierarchical Data Format version 5 (HDF5) files 1 which makes it possible to quickly access the information of the photons that matter for the considered region of the parameter space.
The fitting code is implemented as a user model for the Sherpa fitting package (Freeman et al. 2001;Doe et al. 2007;Refsdal et al. 2009;Brian Refsdal et al. 2011). While the user model uses a Python interface, the work of reading and convolving the raytraced photon and the reflection model data, and interpolating between the simulated grid points is done by a fast C++ code compiled into a Python module with the help of the Boost-Python library (D. Abrahams and R. W. Grosse-Kunstleve) of the Boost distribution 2 .

Fitting Parameters
The KerrC model parameter are listed in Table 2. The parameters include the black hole mass M , the black hole spin a, black hole inclination i, and distance d. The accretion flow is characterized by a single parameter, the mass accretion rateṀ . The model parameter g selects between the wedge-shaped and cone-shaped corona configurations characterized by r C , θ C , T C , and τ C as described above.
The reflection parameters include an amplitude to increase the reflected emission above or below the selfconsistently derived intensity. The other reflection parameters describe the physical conditions in the photosphere of the accretion disk. The parameters include the metallicity A F e relative to solar, the electron density n e at the inner edge of the accretion disk, and a power law index giving how the electron density scales with radial distance (n(r) ∝ r −α ). XILLVER assumes that the photons irradiating the accretion disk have an energy spectrum of disk photons comptonized by electrons of temperature kT e according to the Nthcomp model (Zdziarski et al. 1996;Życki et al. 1999). We consider kT e as a free fitting parameter here. The parameter κ modifies the energy range used for determining the photon index of the emission irradiating the disk. We use κ = 1.5 giving a fitting range from 1.5 keV to 15 keV. As mentioned above, the code calculates all three Stokes parameters. The model parameter χ gives the rotation of the black hole spin axis of the model relative to the celestial north pole (the black hole spin axis of the model is turned counter-clockwise for χ > 0). Finally, the user can exclude photons launched at radial distances outside the l 1 ... l 2 interval from the analysis. This can be used to estimate the impact of disk truncation or shadowing in rough approximation.

Model Evaluation
As described above, we simulated a library for a set of black hole, accretion disk, and corona configurations. When kerrC is called with parameters between the simulated ones, kerrC identifies the nearby simulated configuration nodes, calculates energy spectra for these nodes, and returns an interpolated energy spectrum using linear interpolation in simplices (Weiser & Zarantonello 1988).
For each configuration node, the direct emission data are used to calculate the energy spectrum of the direct emission. Subsequently, the disk data are analyzed to obtain the radial 0.1-1000 keV fluxes f i and the 1.5-15 keV photon indices Γ i of the returning and corona photons irradiating the disk in logarithmicly spaced radial bins (i = 1...12). For the same radial bins, the plasma frame energy spectra F C,i,j (E) of photons scattered into each of 10 inclination bins (j = 1...10) are acquired. The inclination refers to the plasma frame direction of the scattered photon. The subscript C indicates that these energy spectra are based on Chandrasekhar's scattering formalism.
In the last step, the reflected energy spectra are calculated.
For this purpose, each reflected photon enters the analysis with a weight that accounts for the results from the XILLVER radiative transport calculations. We use the energy spectra F X (Γ, A F e , log 10 ξ, kT e , log 10 n e , θ; E) from the XILLVER tables version Cp_3.6 3 . Here, Γ is the photon index of The reflection amplitude should be set to unity to get the self-consistently calculated reflection. c XILLVER requires the photon indices of the emission irradiating the accretion disk. We fit the energy spectra from κE1 to κE2 with κ = 1.5, E1 = 1 keV and E2 = 10 keV. d l1 (l2) give a lower (upper) bound on the radial coordinate of the emission of photons entering the analysis.
the emission irradiating the accretion disk, A F e is the metallicity relative to solar system metallicities, ξ is the ionization parameter, kT e is the electron temperature describing the Comptonized XILLVER input energy spectrum, n e is the electron density of the reflecting plasma, θ is the plasma frame inclination angle of the reflected emission, and E is the plasma frame energy of the incoming and outgoing photon. Our treatment seeks to minimize the impact of the XILLVER input energy spectra by multiplying each photon's statistical weight with the ratio r of the energy spectrum for the actual ionization parameter divided by the energy spectrum for the highest simulated ionization parameter (log 10 ξ max = 4.7): The index i denotes the radial bin r i where the reflection takes place, Γ i and ξ(r i ) are known from the disk analysis, kT e is a fit parameter, n e (r i ) follows from the fit parameters n e,0 and α, and θ and E are the inclination and energy of the reflected photon. The multiplicative correction is justified by the fact that the XILLVER results for (log 10 ξ) max = 4.7 agree with Chandrasekhar's results for a pure electron scattering atmosphere (Garcia 2010, and private communication). The motivation for our treatment is that the reflected Chandrasekhar en-ergy spectrum resembles the reflected energy spectrum to first approximation. The re-weighting factor is used as a correction of the Chandrasekhar result. The treatment leads to more physical results than using the ratio of the XILLVER energy spectrum divided by the incident energy spectrum. In the latter case, the XILLVER assumption of incident power law energy spectra can lead to reflected energy spectra with unphysically high fluxes of high-energy photons, exceeding by far the high-energy photon fluxes provided by the corona. The XILLVER tables are binned in 2999 logarithmicly spaced energy bins from 0.07 to 1000.1 keV (adjacent bins spaced 0.3% apart). We smooth the XILLVER results with a Gaussian kernel with a σ of 6 bins (1.9%) to reduce statistical errors associated with a poor sampling of the XILLVER tables by the simulated photons. Although our reflection treatment conserves the plasma frame photon energy for each individual photon, the reweighing with the XILLVER results accounts for energy gains and losses when averaged over many photons. The code uses linear interpolation in simplices to get the XILLVER energy spectra between the simulated XILLVER nodes. If we hit the edge of the simulated parameter space, we use the edge value.
Our code uses Chandrasekhar's results for the polarization of photons scattering off an indefinitely deep electron atmosphere. For each reflection of a photon, the result accounts for the polarization fraction and direction and the inclination and azimuth angles of the incoming and outgoing photon beams. The polarization treatment is thus rather detailed, but neglects the impact of atomic emission and absorption. As scattering tends to generate polarization, but atomic emission does not, our predicted Fe K-α polarizations are expected to be slightly too high. The treatment furthermore does not account for the difference of the polarization dependence between the Thomson and Klein-Nishina scattering cross sections.
Some XILLVER models did not converge satisfactorily (J. García, private communication). For electron densities n e > 10 19 cm −3 we thus use the XILLVER model for n e = 10 19 cm −3 .
The code saves CPU time by storing and later reusing some of the intermediate results. For example, as long as l 1 and l 2 are not changed, the direct emission results and the disk analysis results can be stored and reused for each configuration node. The reflected emission results need to be re-calculated every time the reflection parameters are changed.

Scaling Relations
The thin disk solution of Page & Thorne (1974) implies the following scaling relations of the temperature scale σ, the observer flux multiplier f i , and the disk flux multiplier f i→j with black hole mass M , mass accretion rateṀ , and source distance d: (16) We extracted these scaling relations by numerically calculating the temperatures and weighting factors for different M andṀ input parameters, and evaluating how they depend on these input parameters. kerrC uses Equation (16) to map the M andṀ values provided by the user to a temperature scale σ, and then uses this σ to interpolate between the simulated σ-values. Equations (17) and (18) are used to adjust the flux multipliers in the analyses of the direct emission, the emission irradiating the disk, and the reflected emission.

Code Validation
We validate KerrC by comparing the predictions for the thermal energy spectrum with the results of KERRBB (Li et al. 2005). Figure 3 shows that the comparison of the kerrC and KERRBB energy spectra typically agree within <∼20%. The most pronounced differences can be recognized at energies below 0.5 keV. The fact that kerrC underpredicts those fluxes can be explained by the fact that it only accounts for the emission from r ≤ 100 r g . At larger radial distances, the accretion disk emits at rather low energies. Some additional differences may stem from KERRBB using a limb brightening proportional to 1 + cos (θ) and kerrC using the limb brightening prescription from Chandrasekhar's indefinitely deep electron scattering atmosphere. Overall, the energy spectra of the two codes agree well.
We have validated the flux energy spectrum and polarization energy spectrum of kerrC with that of the code MONK (Zhang et al. 2019) and found excellent agreement (Zhang et al., private communication).

FIT OF INTERMEDIATE-STATE CYG X-1 DATA WITH KERRC
We show here results from using kerrC to fit the NuS-TAR (Harrison et al. 2013) and Suzaku (Mitsuda et al. 2007) Figure 3. Comparison of kerrC (orange lines) and KERRBB (blue lines) (Li et al. 2005). Panel (a) shows the models for a spin a of 0.98, inclinations i of 27.51 • (left) and 65 • (right), a black hole mass of 21.2 M⊙, an accretion rateṀ of 0.17×10 18 g/sec, a source distance of 2.22 kpc, a hardening factor of 1.8, reflected returning emission switched off, and with limb brightening. The other three panels compare the two models for the i = 27.51 • case of (a) when one of the parameters is changed but all others are held constant, i.e. when the black hole mass is decreased to 15 M⊙ (b), whenṀ is increased to 0.5×10 18 g/sec (c), and when the spin parameter a is reduced to 0.75 (d).
package HEAsoft 6.28. We use furthermore 4,991 sec of Suzaku data aquired with the X-ray Imaging Spectrometer #1 (XIS1 Koyama et al. 2007) on the same days. The NuSTAR and Suzaku data have been published in (Tomsick et al. 2018). The authors kindly shared the Suzaku XIS1 data with us. We only use the 3-20 keV NuSTAR data, as the statistical errors of the kerrC model become appreciable above this energy. Following Tomsick et al. (2018), we use only the 1-1.7 keV and 2.1-8 keV Suzaku XIS1 data.
The estimates of the distance and mass of Cyg X-1 have recently been revised Miller-Jones et al. (2021). The new distance estimate accounts for the impact of an orbital phase dependent attenuation of the radio signal, and localizes the source at a distance of d = (2.22+0.18-0.17) kpc. The revised distance and optical data con-strain the black hole mass to be M = 21.2±2.2 M ⊙ . The orbital plane inclination is 27.51 • with 68% confidence level lower and upper bounds of 26.94 • and 28.28 • .
We fit the kerrC model absorbed with a fixed n H = 0.6 × 10 21 (Tomsick et al. 2018) with the cross sections of Verner et al. (1996) and the elemental abundances of Wilms et al. (2000). We freeze M , d, i, and n H , and assume an accretion disk plasma with solar elemental abundances (A F e = 1).
We fit the kerrC parametersṀ , r C , θ C , kT C , τ C and log 10 n e . We assume that the electron density of the photosphere does not depend on the radial coordinate (α = 0), and we do not investigate the impact of a truncated or shadowed disk. We fitted the data by first evaluating the model at all simulated configuration nodes for one particular set of XILLVER parameters. Subsequently, we started the fitting with Sherpa from the configuration node with the smallest χ 2 -value. We got the smallest χ 2 -values by randomly thawing four parameters at a time and repeating the fit, switching between the levmar and moncar minimization engines. We ran the minimization many times with different starting values. A good number of fits gets stuck in local minima. A wedge-shaped corona gives the best fit with a χ 2 of 7136 for 1582 degrees of freedom (DoF). The best coneshaped corona gives a χ 2 /DoF of 8001/1582. In the following, we discuss the results for the two corona models in turn. The χ 2 /DoF-values are significantly larger than unity showing that systematic rather than statistical errors dominate the error budget. We thus refrain from giving statistical error estimates. A detailed error analysis is outside of the scope of this paper.

The Best-Fit Wedge-Shaped Corona Model
The parameters of the best-fit wedge-shaped corona are listed in Table 3. The model has a modest black hole spin of a = 0.861, lower than the previously published values from the fitting of the thermal state a > 0.9985 (Miller-Jones et al. 2021) and from simi-  larly high values from fitting of the broadened Fe Kα line shapes (Tomsick et al. 2014;Duro et al. 2016;Walton et al. 2016;Tomsick et al. 2018).
The wedge-shaped corona has a half opening angle of 21.8 • and extends from the ISCO at r ISCO = 2.87 r g to 72.9 r g . The corona temperature is 100 keV and the optical depth is 0.256. The photospheric electron density is log 10 n e = 20.8. The XILLVER electron temperature goes to the lowest possible value, indicating that the fit prefers reflection energy spectra predicted for soft energy spectra irradiating the accretion disk. Figure 4 compares the best fit with the Suzaku and NuSTAR data. Overall, the agreement is typically good to a few percent. The modeled 1-1.7 keV energy spectrum is missing a peak in the middle of this energy range. Furthermore, the model does not fully reproduce the Fe K-α peak around 6.4 keV, and is too soft at >15 keV. Figure 5 gives some of the internal kerrC results for the first configuration node chosen by the interpolation engine (a = 0.9, T C = 100 keV, τ C = 0.5, θ C = 45 • , and r C = 100 r g ). This node usually impacts the results most strongly. Panel (a) shows the emission that did not experienced a reflection (blue) and the emission reflected at least once (orange). Both of these components may or may not have experience one or several Compton scatterings in the corona. The reflected emission starts to dominate at and above 8 keV. Around ∼1 keV, the flux without reflection is almost one order of magnitude stronger than the reflected emission. Panel (b) shows the disk frame thermal emission from the disk (dotted orange line), the returning emission (dashed green line),   (15) for the inclination of emission θ = 41.4 • . These energy spectra show the shape of the reflected energy spectra. kerrC chooses for each individual photon the correction factor r according to the photon's reflection angle. and the returning and coronal emission (blue solid line). The dashed green line was calculated for the same accretion disk without a corona. Two interesting conclusions: (i) the thermal flux from the disk approximately equals the coronal flux. Some of the accretion disk photons are energized (Comptonized) in the corona, and come back, giving rise to a coronal energy flux similar to the original accretion disk energy flux; (ii) the photons returning to the disk owing to spacetime curvature alone (without scattering in the corona), make up only a small fraction of the energy flux irradiating the disk. Panel (c) presents the 1.5-15 keV photon indices of the photons irradiating the accretion disk, showing a graduate steepening further away from the black hole. The ionization drops from log 10 ξ ∼2 close to the ISCO to log 10 ξ ∼0 at 100 r g . Finally, the last two panels show the comoving energy spectra of the photons irradiating the disk (Panel (e)) and the emission leaving the disk (Panel (f)) for different radial bins (see the Figure caption for additional details). Interestingly, the fit does not choose a combination of higher n e and lower ξ (which would give more pronounced lines), as the overall shape of the energy spectrum does not fit well for lower ξ-values.

The Best-Fit Cone-Shaped Corona Model
The parameters of the best-fit cone-shaped corona are listed in Table 4. This model has a slightly higher black hole spin of a = 0.921 with an innermost stable circular orbit at r ISCO = 2.17 r g . The corona is rather distant from the black hole and extends from 41.4 r g to 89.3 r g . The opening angle is found at the largest simulated value of 45 • . Overall, the corona is thus far away and rather extended. The optical depth of the corona is with 0.64 larger than in the case of the wedgeshaped corona (0.256), probably to provide the disk with a sufficiently large flux even though the corona is rather distant and lets a good number of photons escape without scattering. The coronal temperature is with 107 keV slightly higher than in the case of the wedge-shaped corona (100 keV). The photospheric electron density is log 10 n e = 16.3, much lower than in the case of the wedge-shaped corona. The XILLVER electron temperature is 50 keV. Similar as for the wedge-shaped corona, the model fits the data within a few percent for most of the range with larger deviations in the 1-1.7 keV energy range, around 6.4 keV, and at the highest energies (Fig. 6). Figure 7 shows detailed information for the first interpolation node (a = 0.95, T C = 250 keV, τ C = 0.75, θ C = 45 • , r 1 = 50 r g , and r 2 = 100 r g ). Panel (a) shows that the reflected energy spectrum starts to dom- Table 4. Best-fit cone-shaped corona model. Only thawed parameters are listed; see Table 3 for frozen parameters (including nH).

Parameter
Result ( inate over the non-scattered disk and coronal emission at much lower energies, i.e. around 2.5 keV. Panel (b) shows that the corona illuminates the accretion disk weakly at r ∼ r ISCO , where the returning thermal emission dominates the disk illumination. The coronal Comptonized emission starts to dominate the disk illumination above ∼5 keV. Panel (c) presents the photon indices of the emission irradiating the disk. The energy spectra soften up to ∼4 r g where the returning emission dominates and harden at larger distances, where the coronal emission dominates. Panel (d) shows the resulting ionization parameter of the disk plasma. The pronounced minimum at ∼ 10 r g results from the combination of the strong evolution of the energy spectrum as a function of the distance from the black hole and the transition from the dominance of the returning radiation to the dominance of the coronal emission. Panel (e) shows how the energy spectra harden with distance from the black hole as the coronal emission gains importance. Panel (f) shows that the reflected energy spectra are almost featureless close to the black hole and only show emission lines for the outermost radial disk bins.
The results indicate that cone-shaped coronae have difficulties Comptonizing a sufficiently large fraction of the accretion disk photons to account for the observed power law emission. The best-fit corona thus resembles a large umbrella-shaped corona with a large opening angle, enabling it to Comptonize a large fraction of the accretion disk photons. The fit chooses furthermore a low electron density and thus a very high disk ioniza- tion which leads to a high-yield of the reflected emission. This result for the stellar mass black hole Cyg X-1 somewhat resembles similar difficulties in some Active Galactic Nuclei (AGNs) where compact coronae close to their black holes cannot explain the observed high reflected fluxes (Dovčiak & Done 2016).

PREDICTED IXPE RESULTS
We developed a code to simulate and fit IXPE observations for any model available in Sherpa, including kerrC. Sherpa and X-Spec user will know the fake command used to combine a model with response matrices to generate simulated energy spectra. Along similar lines, we added a Python object simPol that provides a fake method to generate Stokes Q and Stokes U energy spectra.
The code requires the model of the I, Q, and U energy spectra, the Ancillary Response File (ARFs) (effective detection area as a function of energy), the Response Matrix Files (RMFs) (energy redistribution owing to detection principle and detector effects), and the Modulated Response File (MRF). The MRF is the ARF times the energy-dependent modulation factor µ(E). The modulation factor gives the fractional modulation of the reconstructed polarization directions for a 100% polarized signal and depends on the polarimeter performance, and the event reconstruction methods.
The analysis of X-ray polarization data is based on assigning each detected X-ray photon a Stokes parameter i = 1,q = 2 µ(E) cos (2ψ) and u = 2 µ(E) sin (2ψ) with ψ being the reconstructed polarization direction and E the reconstructed energy of the event (Kislat et al. 2015;Strohmayer 2017). The sums I, Q, and U of the i, q and u values of all events, respectively, are binned in energy and form the basis of the analysis. The data are then analyzed by simultaneously fitting the detected and modeled I, Q and U energy spectra, using the uncertainties √ I on I, and √ 2I/µ on Q and U . The simPol fake method uses the Q and U standard deviations and Q-U covariances to generate observed Q and U energy spectra (Kislat et al. 2015). Polarization fractions p and directions χ can be calculated from p = Q 2 + U 2 /I and ψ = 1 2 arctan2(U, Q). We use the preliminary IXPE ARFs, RMFs, and MRFs from Baldini (2020), and show in the following only the results summed over all three IXPE detectors. Figure 8 presents the polarization energy spectra for the two best-fit models for a simulated 500 ksec IXPE observation of Cyg X-1. The predicted polarization fractions turn out to be rather low for the two models: lower than 1% below 4 keV and about 1% at and above 4 keV. The >4 keV polarization fractions is higher for the wedge-shaped corona than for the cone-shaped corona. In the 4-6 keV energy range the expected polarization directions are parallel to the black hole spin axis. For polarization fractions < 1%, the statistical errors on the predicted Q and U values are not entirely negligible.
The differences between Q/I and U/I energy spectra of the two models come from the different black hole spins, the different polarizations that the photons acquire in the coronae of different shapes, and from the different accretion disk irradiation and reflection patterns. Note that the reflected intensity depends among other factors strongly on ξ(r), which are very different for the two corona models. Figure 9 shows the simulated IXPE results for the two coronae. For both corona models, the detection of a non-vanishing polarization in the IXPE 2-8 keV energy range will be challenging. Note that the actual Q and U values depend on the orientation of the source in the sky. One could consider using a direction reference that aligns the positive Q-axis with the axis of the VLBA jet (Miller-Jones et al. 2021). For the wedgeshaped corona, most of the signal will be expected in Q.

DISCUSSION
In this paper we describe a new X-ray fitting model, kerrC, and its application to intermediate state Suzaku and NuStar observations of the black hole Cyg X-1. We chose an intermediate state observation as the thermal low-energy emission (< 3 keV) constraints the properties of the accretion disk. The power-law and line emission at higher energies (>3 keV) constrains the properties and location of the corona.    Figure 8. Predicted I, Q, and U energy spectra (top), polarization fractions (center), and polarization angles (bottom) for the best-fit wedge-shaped corona model (a) and the best-fit cone-shaped corona model (b). Positive (electric field) polarization angles are measured counter-clockwise from the block hole spin axis. We are looking from above at an accretion disk rotating according to the right hand rule with the black hole spin pointing upwards.  Figure 9. Simulated outcomes of a 500 ksec IXPE observation of Cyg X-1 for the best-fit wedge model (a,c) and the best-fit cone model (b,d). The upper panels show Stokes I and the lower panels show Stokes Q/I and U/I. For Stokes Q/I = 1 and U/I = 0, the electric field of the X-ray emission is parallel to the black hole angular momentum vector; for Stokes Q/I = 0 and U/I = 1, it is 45 • rotated counter-clockwise.
Using kerrC for the analysis of intermediate-state Cyg X-1 observations reveals the following findings: 1. A wedge-shaped corona above and below the accretion disk, and a cone-shaped corona in the funnel regions above and below the black hole describe the shape of the continuum energy spectra adequately, and can produce the observed relative intensities of the thermal and power law emission components.
2. The wedge-shaped corona fits the data better than the cone-shaped corona. We included cone-shaped coronae in the funnel regions above and below the black hole in kerrC as a possible 3-D approximation of the compact lamppost corona model. However, the fit chooses an extended cone-shaped corona with a large opening angle. Even this large cone-shaped corona has difficulties to produce the observed hard emission, and the fit chooses a very high ionization fraction in order to maximize the yield of reflected photons. The high ionization fraction in turn gives a reflected energy spectrum almost entirely devoid of any emission lines.
3. Even the wedge-shaped corona underpredicts the relativistically broadened Fe K-α line. The discrepancies may be reduced by adding model components, e.g., absorption or additional emission caused by material in the system (Tomsick et al. 2014(Tomsick et al. , 2018. We did not explore these options in our analysis. 4. The thermal component constrains the black hole spin, and we obtain spin parameters between 0.861and 0.921, somewhat smaller than the results of Duro et al. (2016); Walton et al. (2016); Tomsick et al. (2018); Miller-Jones et al. (2021). It should be noted that the black hole spin and the location of the inner edge of the disk are somewhat degenerate model parameters. Neglecting the disk truncation can lead to underestimating the spin, as lower spins move the radius of the innermost circular orbit outwards and can thus mimic disk truncations at r > r ISCO . The code neglects the heating of the disk by the returning radiation and by the coronal emission. Including this effect may lower the black hole spin estimate as a hotter disk can mimic a disk extending to a smaller radius.
5. The energy spectra of the emission irradiating the accretion disk are not well described by simple Comptonized energy spectra but have rather complex shapes resulting from the superpositon of returning emission, and partially Comptonized coronal emission (panels (e) in Figures 5 and 7). The result indicates that some of the assumptions underlying inner disk line fitting are too simplistic.
The assumption of power law energy spectrum hitting the disk will be a better approximation in the case of AGNs than in the case of stellar mass black holes, as in the former case lower energy accretion disk photons need a larger number of Compton scatterings before reaching X-ray energies. The larger number of scatterings, and the larger ratio of the energies of the coronal X-ray photons and the disk photons, will both result in a cleaner power law energy spectrum in the X-ray band.
6. For the small inclination of Cyg X-1, the expected polarization signal is small (∼1% at 4-8 keV) but might be detectable with IXPE. These small polarization fractions are a bit smaller than the results from the OSO-8 observations of Cyg X-1 of 2.4%±1.1% at 2.6 keV and 5.3%±2.5% at 5.2 keV (Long et al. 1980). Note that the disk and corona geometries might evolve in time. The state of Cyg X-1 during the OSO-8 observations is poorly constrained.
Possible improvements of kerrC include the generation of a denser grid of simulated configuration nodes and larger numbers of events per node in the parameter regions inhabited by the actual astrophysical systems. It would be desirable to generate and use XILLVER ta-bles for input energy spectra that resemble the simulated energy spectra impinging on the disk more closely.
It would be interesting to estimate the impact of the heating of the disk by the corona, and to improve the treatment of photons returning many times to the disk.
An interesting extension of kerrC would include alternative geometries, including for example hot inner accretion flows in which the inner accretion disk disappears in favor of hot coronal gas, and seed photons enter the hot inner gas from a geometrically thin or thick accretion disk surrounding the hot central region.
It would be desirable to implement the results of polarization dependent radiative transport calculations as in (Taverna et al. 2020(Taverna et al. , 2021Podgorný et al. 2022). However, the Compton scattering cross sections and the scattering matrices for the reflection off the disk are polarization dependent. The approach adopted here of reweighing pre-calculated geodesics would become much more complicated if the plasma properties (e.g. metallicity, ionization parameter) affect the polarization of the emitted and reprocessed photons.
For the near future, our work will focus on using kerrC to fit data from Cyg X-1 and other black holes for a range of different emission states.