The Reliability of Accretion Disk Inclination Derived from X-ray Spectroscopy of Active Galaxies

The inclination angle of substructures in active galaxies gives insights into physical components from scales of the vicinity of the central black hole to the entire host galaxy. We use the self-consistent reflection spectral model \textsc{RELXILL} to measure the inclination of the inner region of accretion disks with broadband ($0.3-78\,\rm keV$) X-ray observations, systematically studying the reliability of this methodology. To test the capability of the model to return statistically consistent results, we analyze multi-epoch, joint XMM-Newton and NuSTAR data of the narrow-line Seyfert~1 galaxy I\,Zwicky\,1 and the broad-line radio galaxy 3C\,382, which exhibit different degrees of spectral complexity and reflection features. As expected, we find that adding more data for analysis narrows the confidence interval and that multi-epoch, joint observations return optimal measurements; however, even single-epoch data can be well-fitted if the reflection component is sufficiently dominant. Mock spectra are used to test the capability of \textsc{RELXILL} to recover input parameters from typical single-epoch, joint observations. We find that inclination is well-recovered at 90\% confidence, with improved constraints at higher reflection fraction and higher inclination. Higher iron abundance and corona temperature tighten the constraints as well, but the effect is not as significant as a higher reflection fraction. The spin, however, have little effect in reflection-based inclination measurements. We conclude that broadband reflection spectroscopy can reliably measure inner accretion disk inclination.


INTRODUCTION
An active galaxy comprises axisymmetric structures on many scales.On the largest dimensions, the principal axes of the host galaxy of the active galactic nucleus (AGN) are anchored by the global angular momentum distribution of the stars, which in general also defines the orientation of the large-scale interstellar medium.At the opposite extreme, the accretion disk around the central black hole (BH) establishes its orientation based on the angular momentum of its latest fueling episode (e.g., Scheuer & Feiler 1996;King et al. 2005;Li et al. 2015).Launched by AGNs, the relativistic radio jet that extends from the vicinity of the central BH to possibly Corresponding author: Rong Du durong@stu.pku.edu.cnbeyond the size of the host galaxy is believed to align with the spin axis of the BH (Blandford & Znajek 1977; however, see Natarajan & Pringle 1998).In the unified model of AGNs, an axisymmetric dusty torus along the equatorial plane of the nucleus obscures the optical or ultraviolet emission from the inner regions (Antonucci 1993).The structure and inclination of the torus can be constrained from detailed studies of the infrared spectral energy distribution (e.g., Zhuang et al. 2018).
The fuelling process of AGNs is expected to affect the alignment of the inner accretion disk with respect to the stellar disk of the host galaxy (e.g., Hopkins et al. 2012).However, in practice, the orientation on the smallest scales is observed to be poorly aligned with the host galaxy on larger scales (e.g., Kinney et al. 2000; see review in Section 3.3.3 of Kormendy & Ho 2013).No preferential alignment exists between AGN jets and the minor axis of elliptical (Birkinshaw & Davies 1985) or Du et al.
disk (Wu et al. 2022) galaxies.The circumnuclear disks traced by water megamasers also appear misaligned with respect to the host galaxy on global scales (e.g., Greenhill et al. 2009;Pjanka et al. 2017).Nonetheless, after discarding a minority of highly misaligned sources, Middleton et al. (2016) report a loose one-to-one correlation, significant at the 3σ level, between the inclination of the galactic stellar disk and the inclination of the AGN accretion disk derived from reflection spectroscopy, hinting on the possibility that AGNs are fed by different mechanisms.
Radiation generated from the accretion disk photoionizes the broad-line region (e.g., Krolik 1998), whose size and structure can be probed by reverberation mapping experiments (Blandford & McKee 1982).Dynamical modeling of the broad-line region reveals that its structure is usually axisymmetric (disk-like; e.g., Pancoast et al. 2014;Li et al. 2018;Williams et al. 2018), although its formation mechanism is still not firmly established.Hypotheses linking accretion disk winds to the origin of the broad-line region have been discussed for years (e.g., Emmering et al. 1992;Murray et al. 1995;Czerny et al. 2017), with recent attention focusing on the scenario of a failed dusty outflow (e.g., Czerny & Hryniewicz 2011;Czerny et al. 2016;Baskin & Laor 2018).Cross-matching the inclination of the accretion disk and the inclination of the broad-line region inferred from dynamical modeling can offer valuable insights into the physical relationship between these two fundamental components of the AGN central engine.
Measuring the inner accretion disk inclination itself plays a crucial role in understanding the environment around supermassive BHs and the processes of accretion and outflow.Theoretically, if we accept general relativity (Einstein 1916), the no-hair conjecture (e.g., Misner et al. 1973, pp. 875-877) limits parameters describing the configuration of every Kerr-Newman BH (Kerr 1963;Newman et al. 1965) to its mass M BH , charge Q1 , and angular momentum J, which can be decomposed further to a dimensionless spin parameter a * := |J|c/GM2 BH and two angular components, including the orientation of the BH.The inclination θ disk of the inner accretion disk in AGNs, defined as the viewing angle of the system with respect to its normal, is believed to align with the orientation of the central BH because of Lense-Thirring precession (the Bardeen-Petterson effect: Bardeen & Petterson 1975;Pringle 1992;Papaloizou & Lin 1995).On a more pragmatic level, knowing the inner disk incli-nation helps to improve estimates of the intrinsic bolometric luminosity of the AGN, assuming that the inner disk aligns with the outer disk emitting in other wavelengths 2 .Disk inclination measurements also help determine the intrinsic values of projected quantities, such as radial velocities.Precise measurements of inclinations would open up the possibility of employing inclination-dependent theoretical models in data analysis.As an instance in which X-ray spectroscopy utilizes the forward-folding technique to minimize the difference between the observed spectra and an instrumentconvolved model, fitting the thermal blackbody continuum of a relativistic, optically thick, geometrically thin accretion disk (Shakura & Sunyaev 1973; extended in the relativistic regime by Novikov & Thorne 1973) with a weak X-ray corona gives loose constraints on BH spin, whose accuracy would be improved with knowledge of the inclination (e.g., Czerny et al. 2011;Done et al. 2013;Reynolds 2021).
Various efforts have been made to measure the inclination angle of AGN accretion disks.For sources that have a jet, an approximate orientation of the beam, which is equal to the inner disk inclination, can be readily estimated from the radio core dominance parameter (e.g., Ghisellini et al. 1993;Wills & Brotherton 1995) via very-long-baseline interferometry measurements of the fastest part of the jet.Other possible approaches include fitting the profile of double-peaked broad Hα emission lines with a relativistic disk model (Eracleous & Halpern 1994;Storchi-Bergmann et al. 1995;e.g., see Storchi-Bergmann et al. 1997for NGC 1097and Ho et al. 2000 for NGC 4450), studying H 2 O megamaser disks (e.g., in NGC 4258; Herrnstein et al. 1999), reproducing nearinfrared and mid-infrared photometric and interferometric observations (e.g., for NGC 1068; Hönig et al. 2007), and measuring the configuration of circumnuclear gas with high-spatial resolution, near-infrared spectroscopy (e.g., for NGC 3227, NGC 4151, and NGC 7469; Hicks & Malkan 2008).In works that assume the unified model of AGNs, the inner disk inclination is deemed to roughly equal the inclination of other substructures, for example, the narrow-line region (Fischer et al. 2013(Fischer et al. , 2014) ) and broad-line region (Wu & Han 2001;Zhang & Wu 2002) clouds.
Another approach to derive key disk properties of AGNs, including inclination, is by fitting the reflection components of their broadband X-ray spectrum (e.g., Fabian et al. 1989).Empirically, the X-ray emission in AGNs shows similar properties, leading to the con-ventional view that electrons from a hot corona in the vicinity of the BH inverse Compton (1923) scatter thermal optical/ultraviolet photons from the accretion disk into the X-rays, creating a power-law continuum with an exponential cutoff (Thorne & Price 1975;Haardt & Maraschi 1991, 1993;Dove et al. 1997;Belmont et al. 2008).Apart from this power-law continuum, observations of AGNs also identify two distinct features commonly depicted as the consequence of coronal X-rays reflected by a standard thin disk (Novikov & Thorne 1973;Shakura & Sunyaev 1973): one is the Fe Kα line at 6 − 7 keV (Fabian et al. 1989;Tanaka et al. 1995), the other the Compton hump peaking at 20 − 40 keV (Guilbert & Rees 1988;Lightman & White 1988;Matt et al. 1991).The former originates from X-ray fluorescence, while the latter traces continuum emission Comptonscattered into the shape of a hump, with the range determined by photoelectric absorption on the red edge and Compton recoil, Klein & Nishina (1929) cross-section, and cutoff temperature on the blue side.As the irradiation process mainly occurs in the inner disk region, gravitational redshift, relativistic beaming, and the Doppler effect modify the emission lines, which stretch the lowenergy wing and create a sharper blueshifted peak.In addition, the relativistically smeared reflection fluorescent lines in the soft band may be a possible origin of the soft excess (e.g., Fabian et al. 2002).On measuring the inclination, the key physical process is the broadening of the Fe Kα profile, which offers information on the velocity of the emitting fluid element as a function of disk inclination and the radius of the emitters.The red wing of the Fe line, used to measure the BH spin (e.g., Reynolds 2014Reynolds , 2021)), is sensitive to the inner radius of the emitters, whereas the blue edge is sensitive to the inclination owing to Doppler boosting (Rees 1966).Consequently, the inclination in disk reflection models is degenerate with the spin and the emissivity profile.The iron abundance of the disk, to the extent that it modifies the strength of the Fe line, also comes into play.Moreover, according to Fabian et al. (2015), broadband X-ray spectroscopy can place stringent constraints on the geometry and radiative mechanisms of the corona.Thus, the measurement of inclination can be influenced by the coronal temperature and the overall proportion of the reflection component (reflection fraction, R f ), which affect the spectrum's flux, the shape of the Compton hump, and the relative strength of the Fe line.
Our larger aim, the subject of forthcoming works in this series, is to explore the connection between the inner accretion disk inclination and other facets of the AGN and its host galaxy.As an initial step, this paper focuses on ascertaining the degree to which the RELXILL reflection model can deliver reliable measurements of the inner disk inclination under conditions that mimic observations of interest.Systematic uncertainties need to be understood given the complexity of the model and the degeneracy of model parameters.Our method is explained in Section 2. Section 3 describes the tests and results on observed spectra, and Section 4 presents the analysis of mock spectra.Section 5 discusses the sensitivity of the inclination measurements to the reflection fraction, iron abundance, coronal temperature, model selection, and some technical concerns about binning.Conclusions appear in Section 6.We adopt a cosmology with H 0 = 70 km s −1 Mpc −1 , Ω Λ = 0.73 and Ω m = 0.27.

Systematic Inclination Measurement
We systematically measure inner accretion disk inclination by analyzing X-ray spectra with RELXILL as a baseline model.Instead of rigidly fitting the spectra by brute-force computation, we exploit the flexibility of RELXILL to evaluate judiciously and systematically models with an increasing number of emission or absorption components (e.g., soft excess and ionized absorbers).The fit jointly solves for θ disk with other parameters, including the spin a * , the iron abundance A Fe , the power-law photon index Γ, the ionization state of the inner disk ξ, the electron temperature of the corona kT e , the reflection fraction R f , and the radial emissivity profile of the primary source (i.e., the corona).We imple-

Du et al.
ment RelxillCp in the RELXILL family to account simultaneously for part of the soft excess, the broad Fe line, and the Compton hump.This choice is motivated by the thermal Comptonization model (NTHCOMP; Zdziarski et al. 1996;Życki et al. 1999) of the corona in RelxillCp, which is more physical than the cutoff power-law model employed in, for example, RELXILL and RelxillD.
To achieve a wide energy range from 0.3 to 78 keV, we combine the high-quality spectra from the pn chargecoupled device (CCD; Strüder et al. 2001) and the Metal Oxide Semiconductor (MOS; Turner et al. 2001) of the European Photon Imaging Camera (EPIC) detectors onboard the X-ray Multi-Mirror Mission (XMM-Newton; Jansen et al. 2001) to capture complex emission and obscuration features in the softer bands, together with spectra from the Nuclear Spectroscopy Telescope Array (NuSTAR; Harrison et al. 2013) mission, which extends the bandpass coverage to include the Compton hump and offers more information on the Fe line.For basic discussions, we always take advantage of XMM-Newton (pn, MOS1, and MOS2) to provide the high-resolution soft X-ray spectrum of the soft excess and complicated absorption components, which are crucial in broadband fitting.
We fit the spectra with RelxillCp in the RELXILL family (version 2.2), as implemented in Xspec (12.12.1;Arnaud 1996).During the fit, we make the underlying assumption that the "disk-related" parameters (θ disk , a * , and A Fe ), except for ionization ξ, do not vary on our timescales of interest, but that the "corona-related" parameters (Γ, kT e , R f , and the emissivity parameters, including the two indices Index 1 and Index 2 , as well as the broken radius R br ), can vary freely.The fixed parameters are tied among the spectra from different epochs; the rest are untied.For simplicity, we lump ξ together with the set of corona-related parameters hereinafter, although it is actually a disk parameter.Our baseline model is a combination of Galactic absorption, a blackbody soft excess, and RelxillCp3 .For different combinations of data, we adopt different technical strategies.
We divide the data and corresponding strategy into five cases: 1. Single-epoch XMM-Newton EPIC.The lack of higher energy observation is the main motivation to use this method.All relevant parameters are set free and fit jointly.
2. Single-epoch XMM-Newton observation and a non-simultaneous NuSTAR observation.We fix θ disk , a * , and A Fe , and fit the data jointly while setting Γ, ξ, kT e , R f , and the emissivity parameters (Index 1 , Index 2 , R br ) free.
3. Single-epoch XMM-Newton observation and a simultaneous NuSTAR observation.In this fairly ideal case, one that is commonly adopted in recent reflection modeling, we tie both the disk-related and corona-related parameters and set free the cross-calibration constant and normalization parameter.5. Multi-epoch, simultaneous XMM-Newton and NuSTAR data.We assume that all parameters in the same epoch are the same but may vary between different epochs.
In all five cases, the blackbody temperature kT is a free parameter but assumed constant within all the spectral epochs, while the Galactic absorption column density is fixed to the value retrieved from the Leiden/Argentine/Bonn Galactic H I Survey (Kalberla et al. 2005).Note that within each simultaneous epoch, the parameters for the spectra from different detectors are tied.Cross-calibration constants (∼ 1) account for the differences between instruments, and normalization parameters are used to compensate for differences in integrated flux and exposure time.Upon achieving the best-fit parameters, we generate a posterior parameter distribution to determine their confidence intervals.

Reliability Tests
In this section, we address the challenges inherent in the RELXILL methodology and outline our design for reliability tests.Despite its self-consistent integration of physical parameters, the model's high dimensionality and nonlinearity pose challenges, particularly in retrieving global minima during spectral fitting.This complexity introduces intricacies in interpreting results, compounded by the multitude of parameter degeneracies involved.
To establish the reliability of our methodology, we systematically discuss three aspects: For a complex model such as this, parameter space exploration and error calculation are non-trivial.We employ a Markov chain Monte Carlo (MCMC) algorithm (Metropolis et al. 1953) to generate the posterior distribution of baseline parameters.We check for chain convergence and use the 90% confidence interval to document the results of spectral fits, facilitating the assessment of the consistency of θ disk measurements across various circumstances.

The Sources
The prototypical narrow-line Seyfert 1 galaxy I Zwicky 1 has a complicated soft X-ray absorption structure and a strong Fe Kα line (Ding et al. 2022), which is thought to be evidence of strong relativistic disk reflection (however, see Reeves & Braito 2019).The warm absorbers (WAs) in the soft X-rays, a combination of fast outflows and corresponding shocked circumnuclear gas (Ding et al. 2022), evolve over long timescales.Changes in ionization may be driven by the different densities of the clumps crossing the observer's line-of-sight, in which the "skin" layer facing the source has a higher ionization state (Costantini et al. 2007;Silva et al. 2018).I Zwicky 1 was observed once in 2020 by NuSTAR in a ∼ 5 day-long exposure and twice by XMM-Newton, each with a ∼ 70 ks exposure (PI: D. R. Wilkins).Because of the strong variability, we divide the NuSTAR data into five epochs according to their overlap with the XMM-Newton observations (as in Ding et al. 2022, see our Table 1).The broad-line radio galaxy 3C 382 exhibits a WA, a soft excess, and an Fe Kα line, but no Compton hump (Torresi et al. 2010;Ballantyne et al. 2014).According to Ursini et al. (2018), the soft excess comes from a warm Comptonization component with a photon index Γ ≈ 2.4 − 2.5, a temperature of 0.6 keV, and a cor-responding optical depth of ∼ 20.These parameters are all in general agreement with those in radio-quiet Seyfert galaxies (e.g., Jin et al. 2012;Petrucci et al. 2018).The Fe Kα line is broad in two epochs and possibly broad in the other three, which, together with the lack of a Compton hump, seemingly undermines the reflection model previously applied to explain the properties of the source (e.g., Sambruna et al. 2011).However, given that the lower limit of the coronal temperature may be as low as 20 − 40 keV and the reflection fraction is merely ∼ 0.1 (Ursini et al. 2018), the reflection model is intrinsically poorly constrained, and additional information is needed to test it.For our purposes, the broadening of the iron line is a sufficient indicator of reflection-based inclination measurements.Compared to I Zwicky 1, 3C 382 has weaker reflection features and less complex outflow structures.It was observed simultaneously five times in 2016 (PI: F. Ursini; Table 3).

Data Reduction
For XMM-Newton data, we extract the light curve and spectrum from the pn and MOS of the EPIC detectors with the System Analysis Software (version 20.0.0) and the calibration files released on 31 March 2022.For each bundle of observational data, we use epproc and emproc to retrieve its event lists, evselect to extract the spectrum of the source and background, and arfgen and rmfgen to generate, respectively, the ancillary files and redistribution matrices.Checking for pile-up and background flares reveals that no corrections are needed for either source.The spectra are extracted from a circular region with a 40 ′′ diameter around the source.In the case of the pn observations, a source-free polygon region     excluding serendipitous sources and chip edges is chosen for background extraction, whereas for MOS observations the background spectrum is extracted from a 300 ′′diameter circular region within the range of the outer CCD.For I Zwicky 1, we combine the spectra from the MOS1 and MOS2 detectors with the task epiccombine to enhance the signal-to-noise ratio (S/N).We rebin the spectra of 3C 382 with specgroup so as to not oversample the full width at half-maximum resolution by more than a factor of 3, and to ensure that each spectral bin contains at least 25 counts.The NuSTAR data are reduced using NuSTAR-DAS version 2.1.2,with event lists retrieved from the two Focal Plane Module (FPM) detectors processed by the standard NUPIPELINE, which is integrated into the HEASoft bundle (6.30.1).We use calibration files from CALDB version 20220802.A source-free 300 ′′ -diameter circle near the source region is used to measure the background.For I Zwicky 1, we extract the spectra and light curve with a small circular region of 30 ′′ diameter centered on the point source, as recommended for faint sources, and then combine the data from FPMA and FPMB following Wilkins et al. (2021).For 3C 382, to stay consistent with XMM-Newton procedures, we also use a circular region of 40 ′′ diameter centered on the point source to extract the spectrum.The source spec-tra of 3C 382 are further rebinned with ftgrouppha in HEASoft such that each spectral bin contains at least 25 counts.

Broadband Model Fitting
Based on the data rebinning scheme, we apply the Cash (1979) statistics (C-stat) for I Zwicky 1 during statistics minimization, and the χ 2 statistics for 3C 382.I Zwicky 1 falls in case 4 and 3C 382 in case 5 discussed in Section 2.1.All the epochs of 3C 382 and epochs b and d of I Zwicky 1 belong to case 3. We produce a case 4 fit for I Zwicky 1 and a case 5 fit for 3C 382.We adopt the conventional scheme by starting from a phenomenological power-law model and gradually modify the model components.Because X-ray spectra are folded, fitting with a simple model first allows us to get an initial overall sketch of the spectra.
We start by fitting the data with a redshifted powerlaw model (zPowerlw) modified by Galactic absorption (TBabs; see Wilms et al. 2000).The broadband spectra of both sources show deviations from the phenomenological model, which require the addition of ionized absorbers at ∼0.5 keV, a broad Fe Kα line, and a high-energy Compton hump, motivating us to use RelxillCp.A mild excess in the residuals of the soft band prompts us to add a redshifted blackbody component (zBBody), but we do not consider intrinsic absorbers as in previous works (Ursini et al. 2018;Ding et al. 2022).As discussed for cases 4 and 5, coronarelated parameters are untied for non-simultaneous observations, while disk-related parameters are tied.Furthermore, we fix parameters of the accretion disk that cannot be derived directly from our analysis, including its inner radius (R in , radius of the innermost stable circular orbit), outer radius (R out = 400 r g , with r g := GM/c 2 ), and density n = 10 15 cm −3 .The broadband model fitting procedure is illustrated in Figure 1 for I Zwicky 1 and in Figure 2 for 3C 3824 .
We remark on a technical issue related to the blackbody component in best-fit models.Figure 1e presents the best-fit model for I Zwicky 1 as constant*TBabs*XSTAR*RelxillCp.In contrast to a direct combination of XSTAR tables and the baseline model (with zBBody), we exclude the blackbody component because the successful modeling of complex      2018) is that our analysis does not include Reflection Grating Spectrometer (RGS; den Herder et al. 2001) data and may not be able to resolve the two states.In terms of fit statistics over degrees of freedom (DOF), the single-state absorber model returns χ 2 /DOF = 5682.6/5004= 1.136, whereas the double-state model yields χ 2 /DOF = 5641.4/4989= 1.131 (∆χ 2 /∆DOF = −41.2/− 15).The improvement is not significant.
We generate the posterior probability distribution for each parameter using MCMC implemented in XSPEC.Adopting the Goodman & Weare (2010) algorithm, we use 200 walkers with a chain of 6 × 10 6 elements (3 × 10 4 iterations).The proposal distribution is generated from a Gaussian distribution around the best-fit values.The chain convergence is tested using graphs to monitor the evolution of parameters and χ 2 with regard to MCMC steps, and with the integrated autocorrelation time τ f (Goodman & Weare 2010) estimated by the function in EMCEE (Foreman- Mackey et al. 2013).We ensure that the MCMC chains are longer than ∼ 1000 τ f , as suggested by Sokal (1996).The first 6 × 10 5 elements are rejected as the burn-in phase prior to storing the chain, in order that the chain forgets the initialization point and samples in the equilibrium state.The results of reflection continuum fitting of I Zwicky 1 and 3C 382 are given in Tables 3 and 4, respectively.The best-fit parameters of the ionized absorbers are shown in Table 5 for I Zwicky 1 (also see Ding et al. 2022, their  Last but not least, Tables 7 and 8 provide the calculated luminosity and fluxes.Cross-calibration constants ∼ 1 yield good agreement between the MOS detectors and PN on XMM-Newton, as well as between the FPMA and FPMB on NuSTAR.

Inclination with Different Combinations
We rearrange the observed data to mimic the different circumstances (cases) discussed in Section 2.1.For I Zwicky 1, the result from epoch b is selected to present a case 3 fit.We produce a case 1 fit with the XMM-Newton observation from epoch b, a case 2 fit with the XMM-Newton spectra from epoch b and NuSTAR spectra from epoch d, and a case 5 fit with epochs b and d.The broadband model fitting in the previous section is a case 4 fit.For 3C 382, the result from epoch c is selected to represent a case 3 fit.We conduct a case 1 fit with the XMM-Newton observation from epoch c, a case 2 fit with the XMM-Newton spectra from epoch c and NuSTAR spectra from epoch b, and a case 4 fit with epochs b and c and NuSTAR data from epochs a, d, and e.The broadband model fitting in the previous section is a case 5 fit.Additionally, we produce case 3 fits for all available simultaneous epochs of both sources.Errors are calculated with MCMC, and chain convergence is assured.We show the inclination measurements of the different data combinations for I Zwicky 1 in Figure 3, while the results for 3C 382 are given in Figure 4.
The inclination measurements of each source are internally consistent, in view of the overlap of the confidence intervals.The consistency of measurements for all the different cases and different epochs confirms the overall reliability of our method.Comparing the results from different cases, it is clear that the more abundant the data, the smaller the errors, which means that more in-formation provided by the broadband spectra does lead to a better fit, hence a more precise measurement.Moreover, if the several groups of spectra are simultaneous, the dimension of the model is reduced, which further diminishes the flexibility and complexity of the model and improves the fits.Examining the results from different epochs reveals that the measured value from the broadband model fitting in Section 3.3 lies within the average value of the five simultaneous epochs.Inclination measurements for both sources can be obtained consistently despite the differences in their outflow structure.Yet, the accuracy of the measurement can be improved if the reflection features in the source spectra are more dominant.For I Zwicky 1 (R f > 1), the inclination measurements from different cases and different epochs only vary a little with confidence intervals mostly overlapping, especially for the two simultaneous epochs, whereas the inclination measurements of 3C 382 (R f ≲ 0.5) fluctuate between different combinations of data.We note that combining more data improves the fits more significantly for 3C 382 than for I Zwicky 1.That said, although case 5 is intuitively the most ideal one among the five cases, the naïve expectation that more data give better performance does not necessarily hold.We have shown that applying RELXILL to single-epoch data that exhibit strong reflection signal (I Zwicky 1) suffices to offer satisfactory results with errors of only several degrees.Yet, for data with weak reflection features like those in the single epochs of 3C 382, the fits can yield dubious results with measured inclination angles varying by tens of degrees.
In summary: inner accretion disk inclinations can be measured more consistently if the reflection component is more dominant in the spectra or if more observations are combined and fitted simultaneously.generated by setting R in to the radius of the innermost stable circular orbit, R out = 400r g , and n = 10 15 cm −3 .The blackbody temperature is set to 0.1 keV, the redshift to z = 0.05557 (van den Bosch et al. 2015), and other parameters are initialized as shown in Table 9.
We study all combinations in Table 9, creating 13,860 distinct simulations.The range for θ disk covers nearly face-on to nearly edge-on systems, which have a high degree of relativistic smearing.The input values of the iron abundance A Fe are solar to super-solar to ensure prominent Fe features and to represent the physical conditions in most observations (García et al. 2013).The photon index Γ and ionization ξ lie within the physical range in AGNs.The thermal temperature kT e of the coronal electrons shapes the high-energy cutoff.Although the physical nature of the corona currently cannot be derived from first principles (for progress in magnetohydrodynamics simulations, see Kinch et al. 2016, 2020and Jiang et al. 2019), a firm upper limit on kT e is given by the lower limit of ∼511 keV imposed by the temperature of electron-positron pair production.Observations typically find coronal temperatures of ∼ 150 keV (e.g., Ricci et al. 2017), and the cutoff energy set in several models is as high as 300 keV to even ∼ 500 keV (e.g., Baloković et al. 2020;Kamraj et al. 2022).
The radial emission profile, which is degenerate with θ disk , is intricately linked to the selection of coronal models-a topic we delve into in Section 5.6.While our primary focus for the mock tests does not involve a comprehensive exploration of the impact of radial emissivity parameters (Index 1 , Index 2 , R br ) on inclination measurements, we assign them reasonable values and conduct tests to evaluate their successful recoveries.
The spin parameter influences the entire simulation grid.We first select a specific value of a * for preliminary discussions on other parameters, and subsequently extend the coverage from a * = −0.998 to +0.998.Theoretically, higher values of a * usually allow for a larger grid for R f and θ disk (see Section 5.4).In the reflection modeling of actual observations, reported spin values often lean toward higher values (a * > 0.9; Reynolds 2021).Therefore, exploring the parameter space with a higher spin is more informative.Moreover, the degeneracies intertwined within a * that affect the retrieval of other parameters are fortunately disentangled with higher spin, as accurate retrieval of the spin is typically achieved for values above 0.8 in mock analyses (Bonson & Gallo 2016;Choudhury et al. 2017;Kammoun et al. 2018).Thus, we initially set a * to its theoretical maximum of 0.998 to mitigate the influence of spin degeneracy and study this more flexible slice of the parameter space.In Section 5.4, we extend the discussion to other values of a *  (Table 9) to assess the effect of spin degeneracy on the precision of inclination measurement.

Mock Analysis
Retrieving results from mock spectra can be challenging.Since we know the input parameters in advance, it is technically equivalent to the circumstance that we have already achieved best-fit values when fitting the observed spectra.Nevertheless, knowing the inputs is an idealized situation, for in practice we have little confi-  the Goodman-Weare MCMC algorithm in XSPEC.Following Choudhury et al. (2017), we set a Gaussian proposal distribution initialized at X + 3σ X , where X is the input model parameter and σ X is the standard deviation estimated from the covariance matrix of a preliminary fit.For boundary inputs θ disk = [3 • , 87 • ] and R f = [0, 10], the covariance information is changed to a diagonal matrix constructed from 100 times the step length of the parameter.For each group of the mock spectra, we use 200 walkers for the chain and sample to a total of 1.2 × 10 6 elements (6000 iterations) with the first 1.2 × 10 5 rejected.As an example of the mock and error calculation, we show a group of spectra in Figure 5 and its corresponding MCMC corner in Figure 6.
We check the parameter recovery of the mock spectra by seeing whether the "measured" (fitted) values are close to the input values, and whether the error range includes the input.Figure 7 compares the measured versus input values of inclination, with colored regions depicting 90% confidence intervals: the dispersion suggests that the input values are generally recoverable.There is little dependence on the iron abundance or coronal temperature.Since our simulation grid is five-dimensional, the figures in the following sections are provided to show the trend of parameter recovery regardless of the specific parameter choice in one or two dimensions.Additionally, larger reflection fraction and inclination both lead to tighter parameter constraints, as to be expected because either large R f or large θ disk increases the relativistic smearing of the intrinsic spectra.Dauser et al. (2014) find a similar trend for the reflection model based on the lamp-post geometry (Matt et al. 1991;Martocchia & Matt 1996;Reynolds & Begelman 1997;Miniutti & Fabian 2004).The confidence interval shrinks most markedly when R f varies from 0.3 to 1, with less dramatic but still notable improvements for R f = 5−10.Not surprisingly, the completely reflectionfree case (R f = 0) can barely measure θ disk ; the posterior distribution (blue region) is very close to a uniform distribution.

Sensitivity to Reflection Fraction
Envisioning the corona as a razor-thin layer that hugs the entire disk, the RelxillCp model defines the reflection fraction as where f disk denotes the fraction of photons hitting the disk and f ∞ the fraction of photons escaping to infinity (Dauser et al. 2016) in the rest-frame of the corona.
Given that the corona precisely overlays the disk with a large covering fraction, we anticipate the reflection fraction to be ∼ 1.Therefore, in our spectral fits, R f is not intended to serve as an indicator of the geometry or relativistic light-bending; instead, it serves as an empirical representation of the relative strength of the reflection component, encompassing the soft excess, the Fe Kα feature, and the Compton hump.Although this parameter lacks physical significance, it is nevertheless important because it provides insight into the proportion of the reflected emission ( This, in turn, helps us to estimate the systematic error and lends confidence in applying the reflection model. Referring back to the mock results, in a typical case of R f = 1 (e.g., Figure 7), the relative error can be as large as 50% to 90% when θ disk decreases to ≲ 60 • .However, for larger θ disk , the error is ≲ 10 • .When R f ≳ 5, the proportion of reflection is already larger than 83%, and further increasing R f does not significantly raise this proportion, resulting in negligible changes to the confidence interval.On the other hand, for R f ≲ 0.3, where the reflection signal is weak (≲ 23%), robust inclination measurements are challenging to achieve.The impact of reflection fraction is also seen in Section 3.3, where the results of θ disk from different simultaneous epochs of I Zwicky 1 are more consistent compared to 3C 382.Indeed, even a single-epoch XMM-Newton ob-   servation of I Zwicky 1 offers a fair constraint on inclination.Without loss of generality, we define a spectrum as reflection-dominant when R f ≳ 5 and reflectionsubordinate when R f ≲ 0.3.In general, it is usually easier to achieve tighter constraints on inner disk inclination for sources with higher R f , and, of course, larger inclination, because the observed spectrum of more edgeon systems naturally has a stronger reflection signal.However, it is also worth noting that combining more groups of data can counterbalance the poor constraints of reflection-subordinate sources, as examplified by the case of 3C 382.

Sensitivity to Iron Abundance
Figure 8 examines the effect of iron abundance on the measurement of inner disk inclination.Large A Fe increases the strength of the whole Fe profile (García et al. 2013).In our simulations, both the absolute value of in-clination offset and the confidence interval of θ disk generally decrease with a larger input value of A Fe .Consistent with the results of previous sections, this trend holds as well for smaller R f , although the constraints on θ disk are significantly looser because reflection spectroscopy is intrinsically less reliable for reflection-subordinate cases.In contrast with R f , the influence of A Fe on the accuracy of spectrum fitting is not as significant, as evidenced by comparing Figure 7 and Figure 8b, or by comparing the effect of R f and A Fe in Figure 8 itself.On the one hand, the inclination results are rather reliable already for sources with solar or mildly super-solar abundance; on the other hand, it implies that the model is not that sensitive to A Fe compared to R f .This, aside from the support for super-solar abundance in AGN reflection modeling (Fabian 2006), also explains why the broadband fitting for 3C 382 returns such a huge value of A Fe with a large error even though we simultaneously fit five epochs.We cannot use A Fe as a basis for judging the robustness of reflection models.

Sensitivity to Coronal Temperature
Figure 9 documents the effect of coronal temperature on the measurement of inner disk inclination.The inclination offset in terms of error generally decreases with larger input values of kT e .Consistent with previous sections, this trend holds regardless of the values of R f and A Fe , although the effect from R f is greater than that from A Fe .Similar to the case of A Fe , the inclinations are fairly reliable for sources with kT e = 30 keV.In all cases, however, |θ disk |/σ < 1, meaning that deviations between the measured and input inclinations all lie within the 1σ confidence.

Spin Degeneracy
Both the spin and the inner disk inclination contribute to the broadening of the Fe Kα emission line, with the red edge of the profile being more sensitive to a * and the blue edge more sensitive to θ disk .Different combinations of a * and θ disk can yield similar Fe profiles, making them challenging to distinguish during spectral fitting.This inherent degeneracy is further complicated by the coupling of both parameters with the entire model.Given the intricate interplay between a * , θ disk , and other model parameters, our mock tests are conducted with specific considerations for two coronarelated parameters (ξ and Γ), disk geometry, and, once again, the physical background of R f : 1.While ξ is better measured with higher values, incorrect measurements of ξ, followed by A Fe , Γ, and θ disk , pose challenges for spin retrieval (Choudhury et al. 2017).We use log ξ = 3 as the input and find that the 90% error of ξ is ≲ 300 erg cm s −1 for all mock spectra.
2. According to Choudhury et al. (2017), underestimating ξ corresponds to overestimating Γ, and vice versa.Lower ξ shifts the spectra toward harder X-rays, as larger Γ produces less soft X-ray emission.Here, ξ influences recombination and photoionization strength, directly shaping the recombination continuum.We generate spectra with ξ = 1000 erg cm s −1 as an input value that can be acceptably retrieved, reducing its influence on a * and the overall model fitting.As the power-law exponent in the spectral model, Γ should be easier to independently constrain in our broadband fits.
3. The value of a * sets a physical lower limit for the inner radius of the accretion disk with a negative correlation, which is the main effect of the spin in observations (Dauser et al. 2010).Besides, changing a * modifies the geodesics in ray tracing, which has a more complicated influence on observed spectra.For a disk geometry, a * = 0.998 supports a disk reaching ∼ 1.24 r g , while a * = −0.998pushes the innermost stable circular orbit to 9 r g .We set the broken radius to 10 r g to avoid physically impossible circumstances of R br < R in .
4. The R in − a * relation also affects the reflection fraction.The maximal possible R f declines for smaller a * or larger R in (Dauser et al. 2014).Although our phenomenological R f is not directly constrained by a * , some extreme combinations in our simulation grid are rare in real observations.5. Aside from the spin, measured R f is further affected by θ disk because the observer only receives flux traveling at a certain orientation with respect to the system.According to Suebsuwong et al. (2006), extremely large R f (≳ 10) is possible only for systems with high inclinations (θ disk ≳ 70 • ).Still, since R f in RelxillCp is not physical, a few extreme combinations are possible in mocks but do not necessarily exist.The only motivation for including those points in the parameter space is to make the discussion more complete computationally.
In conclusion, we provide rationale for selecting the specific values for ξ, Γ, and R br summarized in Table 9.Additionally, our previous analyses indicate that A Fe and kT e have a lesser impact on θ disk derivation compared to R f .Consequently, the ensuing discussions will focus on the interplay between spin and inclination measurements and the influence of the reflection fraction.
0.9 0.5 0.0 0.5 0.9 Input a * 0.9 0.5 0.0 0.5 0.9 Measured a * (a) Input a * 0.9 0.5 0.0 0.5 0.9 Measured a * (b) a * =+0.9 a * =+0.5 a * = 0.0 a * =−0.5 a * =−0.9 0.9 0.5 0.0 0.5 0.9 Input a * 0.0 0.5 Like the case for inclination, the measurement of a * is influenced by its intrinsic value.All three preceding studies on mock spectra of RELXILL focused on the reliability of spin measurements and observed a more accurate recovery of a * with higher R f and larger a * .This general trend is confirmed in our analysis, as depicted Figure 10a.A larger a * allows the inner edge of the disk to approach the BH more closely, resulting in more pronounced redshifting of photons and the formation of a broader Fe line that is easier to detect.Fig-ure 10b examines the effect of inclination on spin measurement, with colored regions depicting 90% confidence intervals.Generally, larger inclinations lead to less constrained spins because the differences in line broadening caused by a * are more evident at smaller θ disk (e.g., see Dauser et al. 2010).The inclination measurements are affected minimally by the spin (Figures 11a and 11b).
Considering the entire simulation grid of 13,860 spectra, 85% exhibit an offset between the input and measured inclination smaller than 10 • , 77% have an incli-nation offset smaller than 5 • , and 62% have an absolute inclination offset smaller than 5 • , with their 90% confidence interval covering the input inclination.All reduced fit statistics are below 1.1.

Comparison with Previous Work
In this section, we compare both observed and simulated spectra analyses with previous work.Overall, eight X-ray studies have been reported for I Zwicky 1 since the advent of XMM-Newton.For 3C 382, four works have been reported with XMM-Newton observations.Three mock analyses have been published for the RELXILL family, all focusing on the recovery of a * .
For I Zwicky 1, Gallo et al. (2004) discovered a broad iron line, a hard X-ray flare, and spectral hardening during the flare with the ∼ 20 ks XMM-Newton observation taken in 2002.Costantini et al. (2007) discovered two WAs in soft X-rays from the RGS data, which was further confirmed by Silva et al. (2018) from their two XMM-Newton observations in 2015.Based on the same observations in our work, Wilkins et al. (2021) and Wilkins et al. (2022) suggested that both the broad Fe Kα line extending to 3 keV and the Compton hump around 25 keV were well explained with the disk reflection model.Using a lamp-post coronal geometry, an ionization gradient, and a twice-broken power-law emissivity, they obtained a corona height of h > 14 r g during the flares and h ≈ 3.7 r g before and after flares.Ding et al. (2022) continued with the idea of reflection and showed that a single broken power-law emissivity implemented by RelxillCp was sufficient, and that an ionization gradient was unnecessary once the absorbers were added correctly (see their Figure 2 and Section 3.2).While most of the reflection parameters (Γ, kT e , R f ) indicated the same evolutionary trend as in Wilkins et al. (2022), their result for the spin (a * > 0.973) is consistent with that (a * > 0.94) of Wilkins et al. (2022); their inclination angle θ disk ≈ 42.3 • is ∼ 10 • smaller; and their R f ≳ 1 is an order of magnitude bigger, which was argued to be the consequence of applying a different raytracing kernel.Modeling the reflection continuum and ionized absorbers as in Ding et al. (2022) yields generally similar results.However, we provide measurements for Index 2 (∼ 3) while they fixed it to 3. Including ultraviolet data from XMM-Newton observations, Rogantini et al. (2022) fit the spectral energy distribution via both a frequentist and a Bayesian method and further acknowledged the components found in previous works.With refl in SPEX (Kaastra et al. 1996) assuming solar metal abundance and zero ionization, they reported a reflection scale factor (equivalent to R f ) of 0.43.Reeves & Braito (2019) analyzed the 2015 XMM-Newton observations (3-10 keV) with four different models consecutively: a P Cygni profile (Done et al. 2007), a photoionization continuum, a reflection model (XILLVER * KDBLUR), and a fast (∼ c), wide-angle disk wind model (Sim et al. 2008).They claimed that the broad Fe K profile and the deep 9 keV absorption can be described by a disk wind, although the reflection component could not be excluded.Their reflection modeling returned an ionization ξ ≈ 10 3 , a reflection fraction ∼ 0.6, a flat single power-law emissivity with an index around 2.2, and an inclination angle even 2-20 degrees larger than that of Wilkins et al. (2022).However, their spectra excluded the hard X-rays, and their model did not have spin as a parameter due to the use of another ray-tracing kernel.On modeling the outflowing wind, they set the wind location to the escape radius of the BH and obtained an accretion disk inclination angle 49.7 • as an indicator of blueshifted absorption features.
For 3C 382, Torresi et al. (2010) resolved a WA in the RGS spectra from the 2008 observation, which was the first observed WA in a broad-line radio galaxy.Sambruna et al. (2011) analyzed the EPIC spectrum from the same observation using parameters derived from spectral fits of joint Suzaku and Swift/BAT data and discovered both a broad and narrow Fe Kα line, a hump above 10 keV, a soft excess, and a WA.They argued for a two-ionization reflection with a highly ionized component reproducing the broad emission line and a mildly ionized component accounting for the narrow emission line, the hump, and the soft excess.Both reflection components contributed to ∼ 10% of the total continuum.Using two XILLVER models, they reported log ξ = 1.54 and 2.93.Due to low S/N, they could not constrain the spin with both RELCONV and RELLINE.However, for an assumed emissivity index of 2.5, the inner and outer radii of the reflection region were around R in ≈ 10 r g and r out ≈ 25 r g , respectively, and the inner accretion disk inclination is around 30 degrees.With the simultaneous NuSTAR and Swift observations in 2012 and the NuSTAR observation in 2013, Ballantyne et al. (2014) found a weak Fe line in the NuSTAR spectra and argued that it originated from the outer regions (≳ 50 r g ) of an outflowing corona and ionized accretion disk.In their attempts to utilize RELXILL, they reported an upper limit for the reflection fraction (R f ≲ 0.1) and a lower limit for the inclination angle (θ disk > 25 • ), which was insensitive to the data.Based on the same observations in our work, Ursini et al. (2018) jointly fit the PN and NuSTAR spectra and argued that the line broadening can be explained by instrumental effects or reflection on the torus or the broad-line re-gion.With RELXILL, they reported a very small reflection fraction consistent with that estimated by Ballantyne et al. (2014), and yet an upper limit of 15 • for the inclination angle.They also calculated the intrinsic line width for all five epochs: σ < 0.2, < 0.35, 0.5 +0.3 −0.2 , 0.13 +0.09 −0.08 , and < 0.2, respectively.Bonson & Gallo (2016) generated over 4000 highquality mock Seyfert 1 spectra to study spin recovery with RELXILL, assuming reflection fractions of 1 and 5.They discovered that most parameters were overestimated with spectral coverage limited to 2.5-10.0keV and that extending the data to 70 keV improved the measurements.However, the model was found insensitive to Index 1 , and a * was well constrained only for large values of a * (e.g., with error ±0.10 for input a * > 0.8).Choudhury et al. (2017) tested the recovery of a * and ξ in the lamp-post flavor of RELXILL-RelxillLpwith 5800 NuSTAR simulations of a bright source.They found that both parameters are well-recovered at 90% confidence with improving constraints at higher R f and a * , and lower source height h.Besides, they stressed the importance of choosing a reasonable initial point for spectral fitting and spanning the parameter space carefully for nonlinear-behaving parameters (e.g., ξ) to avoid the fit being trapped in a local minimum.Kammoun et al. (2018) generated 60 mock spectra with more complex components including warm and neutral absorbers, relativistic (RelxillLp) and distant reflection, and thermal emission.Their blind fit was produced starting with selecting models, which emulated the typical procedure of observed spectral fitting.They reported that neither the absorbers nor the input a * affected the fits significantly, although h was important.Compared with our mock analysis, however, both Bonson & Gallo (2016) and Choudhury et al. (2017) excluded the soft band.Tests in Kammoun et al. (2018), although more complete in spectral coverage, included a plethora of spectral models beyond our main focus.More importantly, all three works generate mock spectra with a long exposure time to show the upper limitations of model performance with high counts and high-S/N data.By contrast, our analysis tests the model performance with data settings more closely related to real observations.In addition, Choudhury et al. (2017) and Kammoun et al. (2018) adopted the lamp-post geometry, whereas we do not assume any particular coronal configuration.

Model Selection
In practice, selecting a model to interpret the X-ray observations of AGNs is challenging because their typically complicated spectra can be fit by several candidate models.The comparisons in the previous section clearly illustrate the statistical and potential systematic uncertainties introduced during model selection (aspect 3 in Section 2.2).For I Zwicky 1, the presence of obvious line broadening and a high-energy hump may seem to provide exclusive evidence for the reflection model, but the wind model cannot be fully excluded.In the case of 3C 382, a weak reflection signal leaves ample room for alternative models, and the results of reflection models are significantly shifted based on data quality and quantity (Section 3.4), as well as the strategy of parameter selection in fitting.Even within the framework of reflection, the choice of the emissivity profile associated with the coronal model also affects the measurement of inclination and other parameters.Therefore, current reflection modeling cannot circumvent case-by-case studies to fit observed spectra.For reflection-subordinate sources, combining multi-epoch broadband spectra and carefully exploring the parameter space becomes crucial, as it can at least minimize the statistical error.
Within the framework of reflection, we choose Relx-illCp over RELXILL or RelxillD due to its thermal Comptonization mechanism of the corona.The main difference among the models in the RELXILL family lies in the configuration of the corona, which might influence inclination measurements because it determines the radial emissivity profile of the system.In practice, the radial emissivity profile can be parameterized for various geometries (e.g., Dauser et al. 2013).For Relx-illCp, the emissivity is given as r −Index1 between R in and R br and r −Index2 between R br and R out .The fact that these empirical parameterizations often find a high Index 1 (≳ 5; see also Section 3.3) has motivated the lamp-post geometry-an isotropically irradiating source located in the vicinity and on the rotational axis of the BH-which has been widely employed in recent years to interpret the X-ray spectra of AGNs (e.g., Wilkins & Fabian 2011;Dauser et al. 2012;Miller et al. 2015;Beuchert et al. 2017).This leads to the advent of the RelxillLp series.We stress, however, that the configuration of the corona remains unclear and difficult to distinguish from spectral analysis (e.g., Tortosa et al. 2018).Further determination of coronal geometry usually requires additional information, for instance from X-ray reverberation mapping (e.g., Reynolds et al. 1999;Fabian et al. 2009;Kara et al. 2016;Caballero-García et al. 2020) or X-ray spectropolarimetry (e.g., Tagliacozzo et al. 2023; but see Ingram et al. 2023).Comparing well-understood, theoretical emissivity profiles to the profile determined from the Fe line, and applying constraints from reverberation lags, can help establish the geometry (Wilkins & Fabian 2012), but this is beyond the scope of this work.Aside from the arguments in coronal geometry, Relx-illLpCp includes another feature compared to Relx-illCp: the ionization gradient.Photoionized by the irradiation from the corona, the accretion disk receives a flux as a decreasing function of the radius, leading to a radial distribution of the ionization parameter, where n is the number density of the plasma receiving an ionizing flux F .In RELXILL, the parameter n (or n H ) usually denotes the hydrogen number density, following XSTAR.However, the n in Equation ( 2) should be the electron number density, equal to 1.2 times the hydrogen density since most H atoms are ionized, and F becomes the net integrated flux in the 1-1000 Ry energy range (F X ; Tarter et al. 1969;García et al. 2013).Thus, the ionization parameter in RelxillCp is ξ = 4πF X /1.2n H . Svoboda et al. (2012) argued that the radial ionization profile should be considered to derive the very steep emissivity profiles (Index 1 ≥ 7), which is degenerate with both the spin and the inclination.
For 3C 382, it is pointless to debate whether we need to include additional features because reflection modeling suffers from a low reflection fraction.But it is worthwhile to see if the ionization gradient affects the inclination angle for typical sources like I Zwicky 1.We conduct model comparison via the likelihood ratio (Λ) and the Akaike (1974) information criterion (AIC).Because we use C-stat for spectral fitting, the likelihoodratio test is equivalent to the χ 2 hypothesis test (Wilks 1938(Wilks , 1963)), where k is the number of model parameters.We apply AIC to compare the ability of models to properly de-scribe the observed spectrum without significant overfitting and under-fitting.The AIC is defined by which not only includes the likelihood function but also punishes the model with more parameters.We calculate AIC for our fit of I Zwicky 1 using RelxillCp as the null model, and we produce alternative model fitting using RelxillLpCp with coronal parameters initialized according to the best-fit values in Wilkins et al. (2022).
During the fit, the density of the disk and the ionization gradient are considered as corona-related parameters, meaning that they are free to vary across spectral epochs.The radial density profile is not considered because its effect is incomparable with that of irradiation.We report the inclination, spin, and other values related to statistical tests in  2022) (R f < 0.5; h > 14 r g to ≳ 3 r g ), even though all parameters evolve in the same trend.
For generality and our purpose of inclination measurements, it is readily acceptable to adopt RelxillCp to use its two emissivity indices (Index 1 and Index 2 ) and the broken radius (R br ), which implicitly has taken into consideration various coronal geometries without further assumptions.However, it is far more difficult to interpret than simply applying a lamp-post geometry with which many of the model parameters are directly calculated (e.g., the relation between R f and θ disk is wellstudied in Dauser et al. 2014).Moreover, the model is insensitive to the inner emissivity Index 1 (Bonson & Gallo 2016), rendering it challenging to be constrained in a blind fit.We suggest assigning proper initial values for the emissivity parameters and always providing the fitted emissivity along with the inclination measurements for reproducibility (as in Section 3.3).

Effect of Soft Excess Interpretation
The interpretation of the soft excess has been a subject of ongoing debate (e.g., Crummy et al. 2006;Done et al. 2012;Jin et al. 2017;García et al. 2019;Gliozzi &  Williams 2020).The soft excess is, at least in part, a reflection feature that influences the determination of parameters such as a * (e.g., Walton et al. 2013), ionization, and the emissivity profile.However, the soft band often introduces complications with complex absorption and scattering phenomena, mimicking the relativistic effects associated with reflection.Consequently, distinguishing whether the soft excess is due to reflection-based relativistic smearing or purely reprocessing poses a challenge.This ambiguity is often addressed using high-S/N broadband spectra (e.g., Guainazzi et al. 2006;Risaliti et al. 2013;Kammoun et al. 2018).In Ding et al. (2022) and our sections on observed spectra, the utility of multi-epoch, broadband, reflection-dominant spectra with well-modeled absorption reduces the complexities.Using 3C 382 as an example, we then delve into a detailed discussion of the phenomenological modeling of the soft excess with little influence of the WA since WAs are well-modeled in Section 3.3.As the Relx-illCp model inadequately describes the soft excess in 3C 382, we choose to model it using an additional phenomenological redshifted blackbody with a temperature tied among spectral epochs (kT = 8.7 +0.3 −0.3 × 10 −2 keV).Despite statistical robustness, this introduces systematic uncertainty into the spin, ionization, and emissivity profile.The blackbody temperature may vary between epochs, and alternatively, adopting the zBBody model may introduce a selection bias.We conduct four spectral fits to elucidate the effect of our zBBody interpretation of the soft excess: I.The spectral fitting in Section 3. II.Untying blackbody temperatures between epochs.III.Ignoring spectral bins with energies lower than 2 keV and directly producing the fit to check for the existence of a selection bias.
IV. Ignoring spectral bins with energies lower than 2 keV and deleting the absorber table components in the soft band.We then directly start an MCMC error calculation with the best-fit parameters from our original fit (Fit I).
Results are presented in Table 11.The inclination θ disk and spin a * show no significant impact in Fit II.Despite introducing four additional model parameters to account for the temperature, the C-stat remains unchanged (∆C-stat = 0.02).The likelihood-ratio test cannot reject the null model (Fit I), and this scenario is slightly disfavored with ∆AIC = +8.The blackbody temperatures across the five epochs align with the tied value in Fit I (|kT a,b,c,d,e −0.087| < 0.002).Fit III yields unattainable results when information in the soft band is completely lost.The values for inclination (θ disk ≤ 3 • ), spin (a * = 0.998), Fe abundance (A Fe = 10), and most emissivity indices (7 out of 10 Indexes( 1,2 ) → 10, 0) are all pegged to their boundaries.The reflection fraction (R f ≲ 0.29) is notably smaller than the values in Table 4. Therefore, addressing the bias of using a specific soft excess model is challenging through a comparison of fits with broadband data and pure hard X-ray data.In Fit IV, we abandon pursuing a fit and instead start error calculation on the known "best-fit" parameters directly.As expected, a * , A Fe , and Γ generally align with Fit I, exhibiting overlapping 90% error.However, R f is smaller (< 0.25), and ionization is larger (log ξ ≳ 3.2).The larger ξ possibly shifts θ disk to a smaller value.Notably, RELXILL still fits the harder energies well in this scenario.
In summary, the epoch-independent blackbody effectively models the soft excess.We emphasize the importance of incorporating information from the soft band to enhance the precision of the reflection model fit.

Data Rebinning
Finally, we address a crucial technical aspect that can introduce biases in spectrum analysis-data rebinning.This process involves redividing and recombining data into a new set of bins.In spectroscopy, the decision to rebin over energy channels occurs after data reduction.While this reduces raw spectrum information, introduces potential errors, and increases the risk of creating artifacts, it is scientifically and statistically justified.
Handling broadband data from diverse instruments necessitates rebinning when assigning equal weights to all instruments.Analyzing spectra with significantly different photon counts, such as combining optical and X-ray spectra or XMM-Newton and NuSTAR spectra with varying X-ray emissions, can be challenging.Highcount spectra would dominate the fit if the raw data were used without proper rebinning5 .Moreover, even a single spectrum from one instrument may lose important information, like emission lines, if data bins are too narrow (e.g., Kaastra 1999;Kaastra & Bleeker 2016).
For statistical convenience, a common practice in Xray observations is to rebin spectra so that each energy bin contains a minimum of 25 counts.This ensures approximately Gaussian data, facilitating interpretation with standard χ 2 statistics.However, blindly adopting this approximation introduces bias in fitting (e.g., Humphrey et al. 2009;Choudhury et al. 2017).Poisson likelihood and Cash statistics are available for unbinned spectra in event-based X-ray observations.Our analysis of mock spectra (Section 4.2) uses Cash statistics.The debate on statistical choices originates from frequentist statistics, and introducing Bayesian approaches to X-ray spectrum analysis can alleviate it.The nested sampling method (Skilling 2004), widely embraced in cosmology (e.g., Planck Collaboration et al. 2020) and gravitational wave studies (e.g., Abbott et al. 2019), could be beneficial (see, e.g., Buchner et al. 2014).This method offers advantages like model selection with the Bayes factor and robust parameter estimation.The nested sampling method is particularly useful for handling multiple modes/maxima in the posterior parameter distribution of reflection models, challenging for typical MCMC algorithms.
Rebinning introduces computational challenges, especially for instruments with high resolution or when the forward-folding process is repeated frequently in spectral fitting with many free parameters (see detailed discussion in Section 2 of Kaastra & Bleeker 2016).This problem is magnified in upcoming high-resolution X-ray missions like XRISM (Tashiro et al. 2018) and Athena (Nandra et al. 2013).While not rebinning the simulated spectra might extend MCMC analysis time, it helps avoid potential biases introduced during data reduction.Given that the simulations only involve a broad Fe line and lack complex emission lines, we believe oversampling does not compromise important information.
In summary, the choice of data binning involves a compromise between scientific and computational con-siderations.Balancing this tension often means rebining spectra to avoid oversampling the instrumental resolution by a factor of 3 (1.5 times the Nyquist rate; see Nyquist 1928;Shannon 1949).The application of systematic approaches like that proposed by Kaastra & Bleeker (2016), considering both S/N and instrumental resolution, needs validation through tests on mock spectra to confirm its ability to reduce computational time while preserving critical spectral properties, especially complex reprocessing features.This remains a topic for future investigation.

CONCLUSIONS
We propose a systematic method to measure the inner accretion disk inclination of AGNs with broadband (0.3 − 78 keV) X-ray reflection spectroscopy.Using RelxillCp from the self-consistent reflection model family RELXILL, without assuming any particular coronal geometry, we test the reliability of measuring inclination by analyzing joint XMM-Newton and NuSTAR spectra of the narrow-line Seyfert 1 galaxy I Zwicky 1 and the broad-line radio galaxy 3C 382, alongside 13,860 simulated spectra that mimic joint XMM-Newton and NuSTAR observations of 3C 382.
Statistically consistent measurement of inclination can be achieved with RelxillCp.More accurate measurements can be obtained if multiple observations are combined and fitted simultaneously, and if the spectrum is more reflection-dominant.Multi-epoch, simultaneous data offer the best constraints on inner accretion disk inclination, but as long as prominent reflection features exist within the relevant instrumental coverage, inclination measurements can still be obtained, even for a single-epoch XMM-Newton observation.The presence of ionizing absorbers, albeit complicating the spectrum analysis process, in principle does not affect the inclination measurements.In the future, a systematic investigation of the impact of warm absorbers on X-ray spectral analysis would be worthwhile.At the current stage, case-by-case treatments of the soft band should be employed to help the fit converge and reduce systematic errors by successful modeling of the absorption.
Our mock tests indicate that systems with higher inclinations and larger reflection fractions recover the input inclinations more reliably.Systems with weaker reflection (R f ≲ 0.3) yield greater scatter in the inclination recovery, whereas for reflection-dominant systems (R f ≳ 5) the precision improves mildly with a larger reflection fraction.The reflection fraction can be treated as a qualitative indicator of the systematic bias of applying RelxillCp.Higher iron abundance and coronal temperature help tighten the constraints on inner ac-Du et al.
cretion disk inclination, but the improvements are not as significant as the reflection fraction.The spin of the BH, although degenerate with inclination, does not affect significantly the measurement of inclination.

Figure 1 .
Figure 1.(a) The XMM-Newton and NuSTAR broadband spectral fit of I Zwicky 1, showing the best fit that contains the Galactic absorbed reflection model, warm absorbers (WAs), and ultra-fast outflows (UFOs).For clarity, only the spectra of epoch d are shown.Note that the spectra of all five epochs are fitted simultaneously.The residuals in terms of sigmas with error bars of size one are shown for various models: (b) the phenomenological redshifted power-law model, (c) RelxillCp, (d) RelxillCp plus a redshifted blackbody component to represent a mild soft excess, and (e) the best-fit model, with ionized absorbers added to the model of panel (c).

Figure 2 .
Figure 2. (a) The XMM-Newton and NuSTAR broadband spectral fit of 3C 382, showing the best fit that contains the baseline model and a warm absorber (WA).For clarity, only the spectra of epoch c are shown.For plotting purposes, the XMM-Newton data have been rebinned to S/N > 30, while the NuSTAR data have been rebinned to S/N > 10.The spectra of five epochs are fitted simultaneously.The residuals in terms of sigmas with error bars of size one are shown for various models: (b) the phenomenological redshifted power-law model, (c) RelxillCp, (d) RelxillCp plus a redshifted blackbody component to represent a mild soft excess, and (e) the best-fit model, with a WA added to the previous model.

Figure 3 .Figure 4 .
Figure 3.A scatter plot showing inclination measurements from (a) different combinations of I Zwicky 1 data groups (cases 1-5) and (b) two simultaneous epochs (b and d).The areas shaded in gray show the confidence range for θ disk derived from broadband model fitting (Section 3.3).The two dashed lines give the average of the upper and lower edges of the 90% error bars of the two epochs.

Figure 5 .
Figure 5. (a) An example of the simulated case 3 spectra, with the model in black lines.(b) The C-stat residuals between the generated spectra and the model.The input values are θ disk = 45 • , a * = 0.998, AFe = 3, R f = 1, and kTe = 150 keV.For plotting purposes, the XMM-Newton data have been rebinned to S/N > 30, and the NuSTAR data have been rebinned to S/N > 10.

Figure 6 .
Figure 6.Results of the MCMC analysis for the relevant parameters of the spectra shown in Figure 5.The red lines correspond to the measured values, and the blue lines give the input values of the simulation.In general, all parameters are well recovered.dencethat best-fit parameters are actual model parameters (global minima), although prior knowledge of the physical properties of the source may serve as a guide.Hence, we need to calculate properly the posterior distribution in mock analysis, such that the sampler explores the parameter space as much as possible.However, for our model with much flexibility and complexity, the sampler can get stuck easily in local minima if we

Figure 7 .
Figure7.Inclination measurements for simulated spectra with solar iron abundance (AFe = 1) and mild coronae (kTe = 30 keV) in the maximum Kerr case (a * = 0.998).The colored regions depict the 90% confidence intervals for different input reflection fractions R f .The dashed line depicts the expected values of θ disk .

Figure 8 .
Figure 8.The dependence of inner disk inclination θ disk on input iron abundance Ae and reflection fraction R f , for fixed kTe = 30 keV, θ disk = 45 • , and a * = 0.998.Panel (a) shows the measurements of θ disk and the 90% confidence intervals, and panel (b) the absolute offset in θ disk .

Figure 9 .
Figure 9.The dependence of the relative offset in inner disk inclination (σ refers to the average statistical error for each point in each panel) on coronal temperature kTe and input (a) reflection fraction R f (with solar AFe) and (b) iron abundance Ae (with R f = 5), for fixed θ disk = 45 • and a * = 0.998.

Figure 10 .
Figure 10.Spin measurements for simulated spectra with different values of (a) reflection fraction R f and fixed θ disk = 45 • , and (b) different values of θ disk and fixed R f = 1.In all panels, we fix AFe = 1 and kTe = 30 keV.

Figure 11 .
Figure 11.Panel (a) shows the dependence of the measurements of θ disk on spin a * with R f = 1; the colored regions depict the 90% confidence intervals and the dashed line the expected values.Panel (b) illustrates the relative offset in θ disk with respect to different a * and R f .In all panels, we fix AFe = 1 and kTe = 30 keV.

Table 1 .
Epochs of I Zwicky 1 Observations

Table 3 .
Best-fit Parameters for the Continuum of I Zwicky 1

Table 4 .
Best-fit Parameters for the Continuum of 3C 382 Epoch

Table 5 .
Best-fit Parameters for Ionized Absorbers of I Zwicky 1

Table 9 .
Parameter Values of the Simulation Grid

Table 10 .
Comparison of Best-fit Results of I Zwicky 1 Using RelxillCp and RelxillLpCp

Table 10 .
(Dauser et al. 2022)nd spin measurements are not significantly affected, the reduced C-stat for RelxillLpCp is slightly larger than for RelxillCp.A likelihood-ratio test cannot reject the null model, and the alternative model is disfavored by ∆AIC = +94.Therefore, the ionization gradient is unnecessary for I Zwicky 1, which is consistent withDing et al. (2022).In addition, returning radiation is also included in RelxillLpCp(Dauser et al. 2022), which is tested jointly.Without a notable effect on θ disk , we have observed stronger reflection (R f > 0.5 to > 4.3) and a lower coronal height (h ≲ 5 r g to ≲ 2 r g ) comparing toWilkins et al. (

Table 11 .
Comparison of Best-fit Results of 3C 382 with Different Treatments of the Soft Band