Inferring the morphology of AGN torus using X-ray spectra: a reliability study

ABSTRACT

Numerous X-ray spectral models have been developed to model emission reprocessed by the torus of an active galactic nucleus (AGN), e.g. UXCLUMPY, CTORUS, and MYTORUS. They span a range of assumed torus geometries and morphologies – some posit smooth gas distributions, and others posit distributions of clouds. It is suspected that given the quality of currently available data, certain model parameters, such as coronal power-law photon index and parameters determining the morphology of the AGN torus, may be poorly constrained due to model degeneracies. In this work, we test the reliability of these models in terms of recovery of parameters and the ability to discern the morphology of the torus using XMM–Newton and NuSTAR spectral data. We perform extensive simulations of X-ray spectra of Compton-thick AGNs under six X-ray spectral models of the torus. We use Bayesian methods to investigate degeneracy between model parameters, distinguish models and determine the dependence of the parameter constraints on the instruments used. For typical exposure times and fluxes for nearby Compton-thick AGN, we find that several parameters across the models used here cannot be well constrained, e.g. the distribution of clouds, the number of clouds in the radial direction, even when the applied model is correct. We also find that Bayesian evidence values can robustly distinguish between a correct and a wrong model only if there is sufficient energy coverage and only if the intrinsic flux of the object is above a particular value determined by the instrument combination and the model considered.

methods: statistical, galaxies: active, X-rays: galaxies

1 INTRODUCTION

The optical spectra of quasars and Seyfert galaxies are subject to broad classification: only Type 1 active galactic nuclei (AGN) display Doppler-broadened lines (Balmer lines, etc.), while both types exhibit narrow emission lines. A gaseous and dusty ‘torus’ was hypothesized to explain this diversity: under an orientation-dependent unification scheme, it blocks the line of sight (LOS) to the central engine and the broad line emission region (Antonucci 1993) in the type 2 objects. For example, one simple depiction of the morphology of the AGN torus was visualized in Urry & Padovani (1995) where it was approximated to be a contiguous axis-symmetric dusty doughnut. In the X-rays, optical type 2s (and a few type 1s) exhibit strong line-of-sight absorption (e.g. Awaki et al. 1991b) which can be attributable to a circumnuclear obscuring torus. This finding implied that the X-ray spectral properties of the AGN were dependent on whether and at which angle the observer’s LOS intersected the torus. Physically such an obscuring torus can be a matter reservoir which feeds the AGN over its duty cycle and additionally can have links to radiatively driven outflows (Hönig et al. 2012). The outflows driving matter away from the central black hole into the galactic environment (Hönig 2019) thus are a potential source of AGN feedback, limiting the mass of the central black hole (Murray, Quataert & Thompson 2005). Thus the different roles the torus plays in the AGN or galactic environment can be correlated with the complexities in its structure. Hence, studying the morphology and nature of the torus is important not just for understanding AGN/Seyfert orientation-dependent unification schemes, but also for luminosity-dependent unification schemes (Ricci et al. 2013), contributions of obscured AGN to the cosmic X-ray background (CXB; Comastri et al. 1995; Risaliti, Maiolino & Salvati 1999; Gilli, Salvati & Hasinger 2001; Treister & Urry 2005; Gilli, Comastri & Hasinger 2007) and potential anisotropic radiative effects on host galaxy processes, such as star formation (Murray et al. 2005; Fabian, Celotti & Erlund 2006) or ionization of diffuse plasma (e.g. Yang, Wilson & Ferruit 2001).

Several studies in IR and X-rays have shed light on the possible structure of the AGN torus. The X-ray study by Awaki et al. (1991a) showed that the iron line Fe K edge could be reproduced sufficiently by a torus model where viewing angle is the only parameter. The sufficiency of a one parameter torus model in Awaki et al. (1991a) could be the result of poor quality of data (in terms of signal-to-noise ratio and degrees of freedom available) analysed. A subsequent study by Alonso-Herrero et al. (2003) discovered that a simplistic torus model with a higher optical depth towards the equator was unable to explain the absence of predicted dichotomy in the steepness in the IR spectra for type 1 and type 2 AGNs. Lutz et al. (2004) and Horst et al. (2006) demonstrated that infrared (IR) emission from dust structures is relatively more isotropic than predicted from a simple continuous-doughnut shape; additionally the ratio of X-ray and mid-infrared (MIR) luminosities is similar for both type-1 and type-2 Seyferts. All these results provide observational evidence against the orientation-only classification of type-1 and type-2 AGN. Also, X-ray eclipse events (e.g. Risaliti, Elvis & Nicastro 2002; Markowitz, Krumpe & Nikutta 2014) produced strong evidence that the obscurer is clumpy, thus challenging the pure orientation-dependent paradigm of AGN classification. Additional evidence for clumpiness in the circumnuclear material comes from MIR spectral fits (Ramos Almeida et al. 2011). It has thus been established that it is not just the orientation but also the nature of the obscurer that determines whether an AGN is seen as type-1 or 2 in optical or obscured or unobscured in X-rays (e.g. Ramos Almeida & Ricci 2017). Elitzur (2008) summarized the observational signatures which indicated the presence of the clumpiness of the torus. An even further complication is a likely dependence of covering fraction on intrinsic luminosity (e.g. Burlon et al. 2011; Ricci et al. 2015).

In this work, our attention is focused towards Compton-thick X-ray obscured AGN and we investigate how accurately the properties of the torus and the central engine can be discerned from their X-ray spectra. The X-ray spectrum of a Compton-thick obscured AGN is comprised of two components viz. the photoelectrically absorbed (zeroth-order continuum) and Compton scattered component (reflection component) of the direct power-law originating in the AGN corona. In addition to the continuum, there are several distinct emission lines from the most abundant elements, notably the 6.4 keV iron line, accompanied by their Compton shoulders (e.g. Ghisellini, Haardt & Matt 1994; Done, Madejski & Smith 1996; Matt 2002; Yaqoob & Murphy 2011). The zeroth-order continuum and the scattered power law has information about the torus morphology and the obscured central engine. Previously in X-ray spectral data fits (e.g. Bianchi et al. 2005), the assumed model for the zeroth-order continuum was an absorbed power law and the model for the scattered component was approximated by using simplistic models viz. reflection from infinite slab (Magdziarz & Zdziarski 1995). Reflection models incorporating a semi-infinite slab are more consistent with an accretion disc rather than a torus. Hence, no physically important parameters corresponding to the torus morphology can be derived when a geometrically inconsistent model (e.g. Murphy & Yaqoob 2009) is used for data analysis. Thus in an attempt to study the physics of the X-ray obscurer, numerous physically motivated models simulating the reprocessed X-ray spectral emission from the AGN torus have been developed over the past decade. These models use complex radiative transfer codes to calculate the X-ray spectrum of a Compton-thick AGN. Some assume a contiguous dust and gas structure; others assume clumpy distribution of dust and gas. Other than the differences in the morphology, the models also assume different aspects of radiative physics like different scattering cross-sections and consideration of some different radiative components. The differences in morphology and/or radiative physics introduce different features in the scattered spectra of the AGN torus e.g. like different nature of the Compton-Reflection Hump (hereafter ‘CRH’), different shape of the iron lines, and Compton shoulders. Most of the simulated models are available as FITS tables and can be used for fitting to real data to estimate the model parameters via spectral analysis software like xspec and isis.

Spectral model fits, including those involving torus models, can be potentially misleading in two different ways. First, some parameters in the given model are degenerate with some other parameters and these degeneracies can affect the quality of simple spectral model fits. From the statistical point of view in complex parameter space, simple |$\rm \chi ^2$|-fit algorithms can get stuck and falsely return wrong values of parameter and uncertainties. Secondly, it might be possible to fit a given data set with different models and obtain good measure of goodness of fit (χ²/dof), but with differing value of fit parameters (Guainazzi et al. 2016) for different models assuming different geometry. Ogawa et al. (2019) have demonstrated similar cases of model degeneracy for un-obscured cases, where both relativistic reflection from the accretion disc, modelled using xillver and the scattered continuum from the AGN torus, modelled using XCLUMPY (Tanimoto et al. 2019) could fit the same data set. This problem of model degeneracy will aggravate as we descend down to low data quality shorter exposures or lower fluxes where the signal-to-noise ratio is poor. The data-model residuals alone can become ineffective in determining the best model for a given data. This brings challenges in measurement of quantities derived from fit parameters; e.g. Eddington ratio (λ_edd) or bolometric luminosity (L_bol) calculated from Γ (e.g. Brightman et al. 2013). So from the viewpoint of data analysis, X-ray astronomers must be informed about which parameters from the models can be constrained correctly, which parameters have the potential to mislead and what are the signatures that can be helpful to realize that the fitting model is not appropriate.

In this paper, we investigate the problem from the viewpoint of an observational astronomer. The goal of our work partially mirrors that of González-Martín et al. (2019) in their testing reliability of AGN IR spectral model fits. We use data simulated using the torus models and the instrument functions of two X-ray satellites, XMM–Newton and NuSTAR, whose combined X-ray coverage spans 0.2–78 keV. We use Bayesian analysis to calculate posteriors and evidence value of a model given a simulated data set(s). The case where data are simulated under and analysed with the same torus model is referred to as the intramodel analysis. The case where data are analysed using a torus model different from the one that has been used to simulate it, is referred to as cross-model analysis. We analyse the behaviour of parameter posteriors and the values of Bayesian evidence for both intra- and cross-model analyses in different conditions to understand their implications on model distinction and parameter determination.

The remainder of this paper is organized as follows: In Section 2, we introduce the spectral models and each of their geometries and parameters, we introduce the Compton-thick model used to test them, and we describe aspects of data simulation and the Bayesian spectral fitting method used. In Sections 3 and 4, we present the results of our model fits, for intra- and cross-model analyses, respectively. In Section 5, we briefly discuss feasibility of detection of a broad discline component in these spectra. In Section 6, we present considerations for exposure time requirements. In Section 7, we discuss our results, including guidelines to the X-ray community regarding making conclusions about the original obscurer morphology and the values of parameter that are estimated from the fits, and caveats associated with our analysis. We summarize our findings in Section 8.

2 METHODOLOGY

In this section, we discuss the models we test, their morphology and the parameters. We also discuss the methodology of our data simulation using the models and our data analysis. This section also contains important abbreviations that we use throughout the work in the subsequent sections.

2.1 Description of models used

In this subsection, we summarize the key features of the different models we tested here.

• MYTORUS: Murphy & Yaqoob (2009) assumes an axi-symmetric doughnut geometry. The model thus implements the classic orientation dependent AGN classification paradigm. The size of the torus is determined by the ratio c/a where c is the distance from the point X-ray source and a is the radius of the torus cross-section. The ratio is c/a = 0.5 (fixed) in the published model. The gas density is uniform throughout the torus. The column density at the equator (N_{H, eq}) is a variable parameter; the LOS column density (N_{H, los}) is derived from the ratio c/a and angle of inclination (θ_i). Here and throughout the paper, we refer θ_i = 0^○ or 90^○ as face-on or edge-on, respectively. The mathematical relation between N_{H, eq} and N_{H, los} is given by

  $$\begin{eqnarray*} N_{\rm H,los} = N_{\rm H,eq}\sqrt{1-\left(\frac{c}{a} \right) ^2 \rm cos^2\theta _{\rm i}} . \end{eqnarray*}$$

  (1)

  The primary radiative processes that reprocesses the input radiation are photoelectric absorption, Compton scattering, and fluorescent line emission. The published model in its most recent version has three separate FITS files: each for the zeroth -order continuum, the scattered continuum and the iron fluorescent lines with their Compton shoulders. The input to the torus is assumed to be a simple power law.

• RXTORUS: This model is based on the radiative transfer code REFLEX (Paltani & Ricci 2017). The geometry assumed is an axi-symmetric doughnut. However, there are important additions to the radiative physics of the model, where photon scattering from electrons bound to the metallic atoms are also taken into account. The scattering cross-sections of bound electrons are modified by binding corrections, which decrease the Compton scattering cross-section and the Rayleigh scattering dominates, mainly at the low energies. This is manifested as excess soft band emission, when compared to the corresponding spectra for MYTORUS. Additionally, c/a is a variable parameter (unlike in MYTORUS). N_{H, los} is a derived parameter calculated from equation (1).

• ETORUS: Ikeda, Awaki & Terashima (2009) assume a continuous torus. This model has a spherical geometry with a biconical cutout at the poles; the cone vertices lie above/below the central point on the symmetry axis. The X-ray source is a point source at the sphere’s center. The ratio of the inner to the outer radius of the torus (r = r_in/r_out) is kept fixed at 0.01 in the published model. The opening angle (θ_o) is variable and for obscured AGN, the angle of inclination (θ_i) has to be greater (θ_i > θ_o). The density of the torus is constant all throughout the volume. The model uses N_{H, eq} as a parameter in data-model fits and the LOS absorption is a function of r, N_{H, eq}, θ_i, and θ_o and is given by

  $$\begin{eqnarray*} N_{\rm H,los} = N_{\rm H,eq}\frac{r \rm (cos \theta _{\rm i} - cos \theta _{\rm o}) + sin(\theta _{\rm i} -\theta _{\rm o}) }{(1-r)(r \rm cos \theta _{\rm i} + sin(\theta _{\rm i} -\theta _{\rm o}))} . \end{eqnarray*}$$

  (2)

  However for our analysis we only test the case where θ_i and θ_o are such that N_{H, los} ≃ N_{H, eq}. The model assumes a Thompson scattering cross-section to calculate the scattered continuum. The published model has only the reflected continuum. The zeroth-order continuum is modelled by a simple zTBABS×CABS×CUTOFFPL in xspec notation, where zTBABS models the direct photoelectric absorption and CABS models the Compton scattering losses. ZGAUSS is used to estimate emission lines.

• BORUS: BORUS (Baloković et al. 2018) assumes a continuous torus with a biconically cutout spherical geometry. Qualitatively, similar to ETORUS but the with the cone vertices coincide with the central point X-ray source. The free parameters and model setup are similar to that of ETORUS. Similar to ETORUS the zeroth order continuum is modelled with a power-law attenuated zTBABS × CABS. BORUS allows for the iron abundance (A_Fe) to be a free parameter in the torus, thus simulating the iron line consistently with the scattered continuum.

• CTORUS: The model (Liu & Li 2014) is based on Geant4 biconical cut-out geometry except the gas distribution is clumpy. The clumps are distributed uniformly between the inner radius R_in and the outer radius R_out, effectively forming a thick shell filled with clouds. The angular limitation of the clump distribution is put by a conical surface for which σ (Liu & Li 2014) is 60^○. σ is not a variable parameter here.

• UXCLUMPY: UXCLUMPY (Buchner et al. 2019) uses the radiative transfer code XARS and cloud distribution used in (Nenkova et al. 2008) given by the mathematical formula:

  $$\begin{eqnarray*} \mathcal {N}(\beta) = \mathcal {N}_0 e^{ -\left(\frac{\beta }{\sigma _{\rm o}} \right)^m} , \end{eqnarray*}$$

  (3)

  where β is the latitude angle. The value of m is set to 2. σ_o sets the width of the torus cloud distribution about the equator. The radial distribution of the cloud is uniform (Nenkova et al. 2008). UXCLUMPY uses such a cloud distribution with an additional optional inner Compton-thick gaseous ring. The Compton-thick ring results in additional absorption in the 7 keV to ∼ 20 keV band, thus increasing the curvature of the CRH. UXCLUMPY uses N_{H, los} as a free parameter and it assumes that the line of sight always intersects with at least one clump irrespective of the inclination (θ_i). However, the distribution of the clouds imply that, as θ_i tends towards the edge or equator (90^○), the number of clumps and hence the probability of obscuration increases.

In Table 1, we summarize the models and their most important features and parameters.

Table 1.

In this table we summarize: the models names, abbreviations that we occasionally use in the parentheses and the corresponding papers (Col. 1), the type of torus and the morphology of gas or clump distribution (Col. 2), energy range set by the upper limit of NuSTAR (Col. 3), their corresponding parameter and the prior range we use in our work for Bayesian analysis (Col. 4), the input values of these parameters we used in our data simulations (Col. 5).

Models	Gas distribution and morphology	Energy range (keV)	Free parameters and ranges	Parameter input
(1)	(2)	(3)	(4)	(5)
MYTORUS (MYT)	Continuous:	0.6–78	|$N_{\rm H,eq}(\rm 10^{22} cm^{-2}) = [1,10^3]$|	150(MCT), 680(HCT)
(Murphy & Yaqoob 2009)	Classic Doughnut		Γ = [1.6, 2.5]	1.9
			θi(°) = [0, 90](θi < 60°)	70
			Csc-pl = [10−7, 1]	10−3
			T/R = [0.5, 5.0]	1, 1.8
RXTORUS (RXT)	Continuous:	0.6–78	|$N_{\rm H,eq}(\rm 10^{22} cm^{-2}) = [1,10^3]$|	150
(Paltani & Ricci 2017)	Classic Doughnut		Γ = [1.6, 2.5]	1.9
			θi(°) = [0, 90](θi < 60°)	70
			c/a = [0.1, 1]	0.5
			Csc-pl = [10−7, 1]	10−3
			T/R = [0.5, 5.0]	1, 1.8
ETORUS (ETOR)	Continuous:	1–78	|$N_{\rm H,eq} \simeq N_{\rm H,los} (10^{22}\rm cm^{-2}) = [1,10^3]$|	150, 500
(Ikeda et al. 2009)	Sphere with biconical cutout		Γ = [1.6, 2.5]	1.9
	cone vertex above/below the		θo(°) = [0, 70]	45
	axis X-ray source on the symmetry		θi(°) = [0, 90](θi > θo)	60
			Csc-pl = [10−7, 1]	1E − 3
			T/R = [0.5, 5.0]	1
BORUS (BOR)	Continuous:	1-78	|$N_{\rm H,los}(\rm 10^{22} cm^{-2}) = [1,10^3]$|	100(MCT), 500(HCT)
(Baloković et al. 2018)	Sphere with biconical cutout		Γ = [1.6, 2.5]	1.9
	cone vertex co-incident		|$E_{\rm cut}(\rm keV) = [60,500]$|	400
	with the X-ray source		|$\rm C_{\rm frac,tor} = cos \theta _{\rm o} = [0,1]$|	0.58
			|$\rm cos \theta _{\rm i} = [0,1]$| (⁠|$\rm cos \theta _{i}\lt \rm cos \theta _{\rm o}$|⁠)	0.48
			AFe = [0.1, 10]	1.0
			Csc-pl = [10−7, 1]	1e−3
			T/R = [0.5, 5.0]	1
CTORUS (CTOR)	Clumpy:	1.2-78	|$N_{\rm H,los}(\rm 10^{22} cm^{-2}) = [1,10^3]$|	100(MCT), 500(HCT)
(Liu & Li 2014)	Uniformly distributed		Γ = [1.4, 2.5]	1.9
	clouds in a thick spherical		Ncloud = [2, 10]	4
	shell		cosθi = [0.05, 0.95](θi > 30°)	0.35
			Csc-pl = [10−7, 1]	10−3
			T/R = [0.5, 5.0]	1,1.8
UXCLUMPY (UXCL)	Clumpy:	0.3-78	|$N_{\rm H,los}(\rm 10^{22} cm^{-2}) = [1,10^4]$|	100(MCT), 500(HCT)
(Buchner et al. 2019)			Γ = [1.4, 2.5]	1.9
	Distribution proposed in		|$E_{\rm cut}(\rm keV) = [60,500]$|	200,400
	(Nenkova et al. 2008)		σo(°) = [6, 84]	45
			Cfrac = [0.0, 0.6](inner ring)	0.4
			Csc-pl = [10−7, 1]	10−2(MCT)
				3.5E − 3(HCT)
			θi(°) = [0, 90]	60

2.2 The implemented model

We simulated the data while implementing the same basic model components across all the torus models. For simplicity, we choose one representative model, typical for X-ray obscured Seyferts in the nearby Universe. The model contains the following components:

• Zeroth-order continuum: The intrisic spectrum of the corona is assumed to be a simple power law or a cutoff power law (ICPL hereafter). In case of an obscured AGNs, the ICPL from the AGN corona is attenuated by the column of the AGN torus, by the process of photoelectric absorption. This absorbed ICPL is referred to as the zeroth-order continuum (Murphy & Yaqoob 2009) (I_{tor, t}(E, N_{H, tor}, Γ, θ_i, R, x₁, x₂,..), where x₁, x₂,.. are morphological parameters of the torus). The parameters that directly affect the zeroth-order continuum involve the parameters of the ICPL viz. Γ, E_cut and the N_{H, los} column of the torus. Throughout the paper N_{H, los} and N_{H, eq} are expressed in units of 10²² cm⁻² unless stated otherwise. For the models MYTORUS, RXTORUS, CTORUS, and UXCLUMPY, the zeroth-order continuum is published as FITS tables, whereas models like BORUS and ETORUS do not have a published FITS table for the same and hence one uses ZTBABS × CABS¹× ZPOWERLAW to simulate it.

• Scattered/reflected continuum: In addition to photoelectric absorption, the incident X-ray photons undergoes Compton scattering (once or multiple times) in the torus material. This modifies the ICPL to give rise to the scattered continuum, I_{tor, R}(E, N_{H, tor}, Γ, θ_i, R, x₁, x₂,..), with distinct features like the CRH, fluorescent emission lines and a softer tail (E < 6 keV) which are discussed below. For all the models the scattered component are published FITS tables. The scattered continuum contains the information about torus morphology. One aspect of the application of a torus model is the value relative normalization (denoted as T/R hereafter) of the zeroth-order continuum with respect to the scattered continuum. The user has the freedom to either set T/R = 1 or let T/R vary. In a real scenario the scattered component is delayed with respect to the zeroth-order continuum due to light traveltime. Thus, in such a case if the ICPL exhibits significant variability T/R ≠ 1 (see MYTORUS online manual²) and should be kept free in the analysis.

• Emission lines: The spectra also contains several fluorescent emission lines, of which the most notable is the iron emission line complex. The lines are accompanied by their corresponding Compton shoulders, which is the resultant of higher order down-scatterings. The treatment of lines are different in different models, e.g. MYTORUS provides only two emission lines, whereas UXCLUMPY, CTORUS, RXTORUS, BORUS provides numerous emission lines. ETORUS does not provide any lines, so the user needs to add an external model e.g. zgauss to simulate them.

• Scattered power-law/warm mirror reflection: Many heavily absorbed AGNs show a significant excess emission component in the soft band. This component might be the scattered component from diffuse gas which cannot be obscured by the torus (e.g. Bianchi, Guainazzi & Chiaberge 2006; Brightman et al. 2014; Buchner et al. 2014, 2019) or the scattered component from the volume filling interclump medium of a clumpy torus (Buchner et al. 2019). For simulation of data using MYTORUS, RXTORUS, BORUS, ETORUS, and CTORUS the scattered power law or warm mirror (both referred as SCPL hereafter) was assumed to be a simple or cutoff power law. UXCLUMPY allows two setups for data simulation and fitting. In one setup, the user can use a simple zcutoffpl for the SCPL. We call this setup UX-ZCPL hereafter in the paper. In another setup, where the SCPL is the scattered emission from the interclump medium of a clumpy torus where scatterings of several orders are considered (Buchner et al. 2019), the user can use the published UXCLUMPY model component, uxclumpy-cutoff-omni. As a result, this component contributes its own weak reflection hump and Fe Kα line. We call this setup UX-OMNI hereafter in this paper. The warm mirror spectrum thus emits its own mild CRH and a weak iron line. The ratio of the normalization (C_sc-pl) of the SCPL with respect to the scattered component ranges from 10⁻³ to 10⁻² times that of the scattered component for all simulations. In this work, we refer to the ratio of the scattered power-law normalization to the torus normalization as C_sc-pl.

• Soft X-ray thermal emission: We can expect that X-ray spectra will almost always contain some potential contamination from host galaxy structures such as star-forming regions and point sources such as ULXs or XRBs, or cluster gas if one is studying an AGN in a cluster or X-ray-emitting shocks for radio-loud objects. To encorporate such soft contamination by significant emission lines caused by photonization (Kinkhabwala et al. 2002), we introduce two components of ionized gas emission. Correct modelling of the component would involve use of e.g. photemis³ from the package xtardb. This is however computationally prohibitive for the simulations we produce. Instead, we use collisionally ionized plasmas (apec;(AtomDB version 3.0.9 ⁴). The different shape has no impact on how torus models are fit.

• Galactic absorption: We use TBABS (Wilms, Allen & McCray 2000) to simulate the Galactic absorption. We assume a value of N_{H, Gal} of 10²¹ cm⁻² for all the simulated data.

The generic model of Seyfert-2 (see Fig. 1) implemented for all the available torus geometries and morphologies is

$$\begin{eqnarray*} f(E)&= & e^{-\tau (N_{\rm H,Gal},E)}[f_{\rm apec}(E,T_1)+f_{\rm apec}(E,T_2)+I_{\rm sc-pl}(E,\Gamma) \nonumber \\ && + I_{\rm tor,T}(E,N_{\rm H,tor},\Gamma ,\theta _{\rm i}, T ,x_1,x_2,..) \nonumber \\ && + I_{\rm tor,R}(E,N_{\rm H,tor},\Gamma ,\theta _{\rm i}, R , x_1,x_2,..) \nonumber \\ && + I_{\rm tor,lines}(E,N_{\rm H,tor},\Gamma ,\theta _{\rm i}, R , x_1,x_2,..)] . \end{eqnarray*}$$

(4)

The regime of Compton-thick AGNs practically starts from N_{H, los} ≃ 100 where the scattered continuum contributes a significant amount of flux to the spectra in addition to the zeroth-order continuum. If we move to higher values of the LOS column density, the zeroth-order continuum drops drastically and the scattered continuum starts dominating. At regime where N_{H, los} < 200 the zeroth-order continuum is dominant in the E > 4 keV band compared to the scattered continuum. In the regime of N_{H, los} > 200 the scattered continuum starts overwhelming the zeroth-order continuum. So in context of this work, we phrase the regime with N_{H, los} < =200 as medium Compton-thick (MCT) and N_{H, los} > 200 as heavy Compton-thick (HCT) regime. In the heavy Compton-thick regime we analyse two different classes of spectra for two models (CTORUS and BORUS). The classes are referred to as HCT0 and HCT1. The scattered power law (SCPL) is comparatively stronger in the HCT0 class compared to the HCT1 class in the 2–10 keV energy range (see Fig. 2d and f). Its contribution is weaker with respect to the total flux (⁠|$F_{\rm 2-10,torus}/F_{\rm 2-10,total} \le 70{{\ \rm per\ cent}}$|⁠) in the HCT0 case and the same ratio being stronger (⁠|$\gt 70{{\ \rm per\ cent}}$|⁠) in HCT1 case for each of two models, CTORUS and BORUS.

Figure 1.

The plot illustrates the model components and how they contribute to the total spectrum. The plot is generated under UXCLUMPY with N_{H, los} = 100, which implies a medium Compton-thick (MCT) regime, and it is clear that the zeroth-order continuum component has a higher flux contribution in the E > 5 keV band.

Figure 2.

We plot overplot the spectrum of the input models for both the MCT and the HCT regime. HCT0 refers to the spectra where the CRH is comparatively weaker compared to the HCT1. For HCT0, the quantity |$r = F_{\rm 2-10,torus}/F_{\rm 2-10,total} \lt 70{{\ \rm per\ cent}}$| whereas, for HCT1 it is |$\lt 70{{\ \rm per\ cent}}$|⁠. In our cases CTORUS and BORUS, |$r_{\rm HCT1} \simeq 90 {{\ \rm per\ cent}}$| for both and r_HCT0 ≃ 65 per cent and 45 per cent for BORUS and CTORUS, respectively.

2.3 Data simulation

Data were simulated with the fakeit command in xspec for XMM–Newton EPIC-pn (Jansen et al. 2001) and NuSTAR focus plane modules (FPM) A and B (Harrison et al. 2013) based on the instrument responses. For our primary analysis, we assume joint simultaneous observations with both missions, simulating good exposure times after screening of 100 ks for XMM–Newton (t_XMM) and 50 ks per NuSTAR FPM module (t_NuSTAR); later, we explore the effects of using only one mission on the model fits. In our primary analysis, the 2–10 keV flux (f_2 − 10) of observed and absorbed flux for all models/instruments was ∼0.5 mCrb. However, if specific cases require us to analyse a spectrum where 2–10 keV flux is significantly different from 0.5 mCrab i.e. |f − 0.5| ≥ 0.09 we specify the values of flux explicitly. For our primary analysis, we simulated data for 100 ks (t_XMM) and 50 ks (t_NuSTAR) on the XMM–Newton EPIC-pn and NuSTAR instruments, respectively, when the analysis concerned the behaviour of the parameters. In a separate study (the results of which are presented in Section 3.3), we intend to understand the distribution of statistical properties (e.g. the 90 per cent confidence region) of the parameters posteriors. To study the distribution of posterior properties we simulate 100 different spectral data sets under a select model for independent fitting using BXA.

For the purpose of studying the dependence of the Bayesian evidence (Z) on intrinsic flux of the object (see Section 6 below), we simulated 10 spectra with the 2–10 keV flux set to be f = 0.83/n mCrb, where n runs from 1 to 10 with the exposure fixed at t_XMM = 100 ks and t_NuSTAR = 50 ks.

2.4 Fitting methods

We carry out Bayesian analysis on the simulated data sets for a given torus model to calculate the parameter posterior distribution. In this paper, we will frequently use common terminology used in simple (e.g. least squares) data-model fitting, but it should be clear to the reader, that despite using the terms that may also apply to simple data-model fitting, the process is Bayesian analysis unless explicitly mentioned otherwise. Specifically, we use nested sampling using multinest (Skilling 2004; Feroz, Hobson & Bridges 2009) implemented via BXA and PyMultinest packages (Buchner et al. 2014) for X-spec version 12.10.1f (Arnaud 1996) for posterior calculation. The issues related to convergence e.g. getting trapped in a local maximum of the likelihood, in standard Goodman–Weare Markov Chain Monte Carlo (GW-MCMC) calculations are not present in nested sampling algorithms. Another convergence problem we have with MCMC is assumption on the length of the chain, which might not be enough for convergence and might require multiple burn-ins and several sequential runs. Nested sampling algorithms, including multinest, attempt to map out all of the most probable regions of parameter sub-space: it maintains a set of parameter vectors of fixed length, and removes the least-likely point, replacing it with a point with a higher likelihood, and thus shrinking the volume of parameter space in each calculation. The convergence of a multinest run is automatic and does not require any initial assumption like the chain length or burn-in. multinest calculates a parameter called evidence Z which is the probability of data (D) given a hypothesis (H), P(D|H), mathematically expressed as

$$\begin{eqnarray*} Z = \int \mathcal {L}(D|\Theta ,H) \pi (\Theta |H) d^N \Theta . \end{eqnarray*}$$

(5)

From a probalistic point of view the evidence (Z) is the probalility of data given the hypothesis, P(D|H). There are two contributors to the final value of evidence (Z) or P(D|H) namely, the likelihood function (⁠|$\mathcal {L}$|⁠) and the final posterior (π) (see equation 5). The likelihood probability here quantifies the deviation of the data from the model and the posterior distribution π(Θ|H) quantifies the volume of the parameter space which signifies the effective implementation of Occam’s razor criterion. The evidence value is a combined effect of these data-model deviations and Occam’s razor. We use multinest version 3.10 with default arguments (400 live points, sampling efficiency of 0.8) set in BXA version 3.31. For all model fits we use uniform or log-uniform as our initial prior parameter distribution. We use cstat (Cash 1979) in xspec as the likelihood function for all model fits. The cstat likelihood is mathematically expressed as |$\mathcal {L}(\Theta ,D|H) = \prod _{i} \frac{(f_i(\Theta) t_{\rm exp})^{D_i}}{D_i!}e^{-f_i(\Theta) t_{\rm exp}}$|⁠, where t_exp is the exposure time for a given instrument, f_i(Θ) is the theoretical count rate for a model H folded with instrument response, and D_i the count rate at the ith channel or energy bin.

We consider the best-fitting value of the parameters to be the median or 0.5th quantile of the posterior distribution and the lower and upper bound of errors are the 0.05th and 0.95th quantile of the distribution, unless stated otherwise. Our prescription of best-fitting value and bound on the error works best when th posterior distribution is a monomodal distribution. Interpretation of the results for the posteriors which show multiple maxima i.e. having multiple solutions, depend on the specific situation concerning the model and the data set. Posterior distributions which are strongly non-Gaussian (i.e. has a minima, are uniform or those that converge towards the edge of the prior range or any other ‘irregularities’) will be termed in general as irregular posteriors hereafter. Interpretation of irregular posteriors will be contextual. Our prescription of best-fitting value and the errors might will be less reliable in such cases and a separate treatment might be required. To quantitatively express the goodness of the recovery of the parameters, we define some simple mathematical parameters to quantify constraints, discrepancy, and the quantitative estimate on the accuracy of the returned values of the parameters.

$$\begin{eqnarray*} R_{90} = q_{0.95} - q_{0.05} \end{eqnarray*}$$

(6)

$$\begin{eqnarray*} \delta = \frac{0.5 R_{90}}{q_{0.5}} \end{eqnarray*}$$

(7)

$$\begin{eqnarray*} \Delta q = |q_{0.5}-q_{\rm input}| \end{eqnarray*}$$

(8)

$$\begin{eqnarray*} r_{\rm q} = \frac{2 \Delta q}{R_{90}} . \end{eqnarray*}$$

(9)

Here q_r is the r-th quantile of the posterior⁵ distribution and thus R₉₀ quantifies the spread of the posteriors in the 90 per cent confidence region. δ is the proportional error about the median, thus the lesser the value of δ the lesser the spread of the 90 per cent confidence region about the median. Δq quantifies the deviation of the median value of the posterior distribution from the input. r is the ratio of the deviation (Δ) to the 90 per cent confidence range (R₉₀). r = 0 indicates no discrepancy, whereas r ∼ 1 indicates that the discrepancy is almost equal to the average statistical error. If r > 1 then the discrepancy is more than the average statistical error. In this work, we take 1.5 to be the critical value for r, thus we term a parameter recovered if 0 ≤ r ≤ 1.5.

3 RESULTS: INTRAMODEL FITS

In this section, we present the results for the intramodel fits (IM-fits hereafter) for joint XMM–Newton EPIC-pn and NuSTAR FPM-A and B data. We grouped all simulated spectra to 30 cts bin⁻¹ for our primary analysis.

3.1 Medium Compton-thick regime (MCT)

In this subsection, we discuss the results of intramodel fits to the simulated data in the MCT regime. Some models adopt LOS column density (N_{H, los}) like BORUS, UXCLMPY, and CTORUS. Other models adopt column density at equator (N_{H, eq}) like MYTORUS, RXTORUS, and ETORUS. During simulation of the data, for the models which adopt N_{H, los} as a parameter, we have set N_{H, los} = 100. For the models that use N_{H, eq} as a parameter, we have set N_{H, eq} = 150 where the implied value of N_{H, los} will be lower than 150, for θ_i situations where the LOS interferes with the torus dust structure (e.g. θ_i > θ_o).

• Column density (N_H): Both LOS column delsity (N_{H, los}) in the models BORUS, CTORUS, and UXCLUMPY and equatorial (N_{H, eq}) column density in MYTORUS, RXTORUS, and ETORUS are recovered with very tight constraints for all models. Over all the tested cases and models, N_{H, los} is recovered, with maximum value of δ = 0.032 and the average values of δ = 0.023. In this regime, the zeroth-order continuum contributes to the total observed flux more significantly than the CRH (⁠|$\left[\frac{F_{\rm T}}{F_{\rm R}} \right]_{\rm 3-100 keV} \sim$| 1.8–5.0). The Fe K absorption edge and the rollover at E < 5 keV are the distinct features of a zeroth-order continuum determining N_{H, los}. The geometry of the system does not play a very significant role in determining the N_{H, los} in the MCT regime, both in terms of tightness of the posteriors and parameter recovery. This holds in the cases where the spectrum is dominated by the zeroth-order continuum.

• Photon index (Γ): We have simulated the data with SCPL-Γ (scattered power law) and ICPL-Γ (coronal power law) set to the same value (Γ = 1.9) for all input models. When fit with the SCPL and ICPL Γ tied, Γ was recovered for all models with values of R₉₀, δ, and r_q never exceeding 0.12, 0.03, and 0.92, respectively. Despite the degeneracy or correlation between Γ and E_cut for (BORUS and UXCLUMPY), the constraints on Γ are tight. We also test the case where the values of SCPL- and ICPL-Γ were not tied for three models. The trends were model dependent in this case. For models CTORUS and MYTORUS (no E_cut) constraints on ICPL-Γ were unchanged. However, the behaviour of SCPL-Γ worsened in both the cases as posteriors widened (R₉₀ = 0.27 for CTORUS, R₉₀ = 0.11 for MYTORUS). In UXCLUMPY for data simulated in the UX-OMNI setup, R₉₀ value of ICPL-Γ worsened for both the UX-OMNI and UX-ZCPL setup due absence of leverage from the SCPL and a simultaneous varying E_cut, in the uncoupled configuration. The trends in SCPL-Γ remained unchanged in the UX-OMNI setup. In the UX-ZCPL setup, the recovered value of SCPL-Γ was flatter (1.83 ± 0.05) and was just barely recovered (R₉₀ = 0.1, r_q = 1.4). For this case, the SCPL-Γ lowered to account for the excess flux (noticeable in the 2 < E < 5 keV band) in the data (simulated in the UX-OMNI setup). Thus it can be concluded that the behaviour of SCPL-Γ can affect constraints on ICPL-Γ when they are tied, more so when they have dissimilar values.

• High-energy cutoff (E_cut): For most of the cases, we simulated the data sets with E_cut set to the highest possible limit allowed by the models, to have a minimum influence of E_cut on the other features (e.g. CRH) in the E < 78 keV band. When we keep E_cut free, the posteriors of E_cut are broad and the constraints are wide in both BORUS and UXCLUMPY. We also simulated two data sets in BORUS and UXCLUMPY with E_{cut, input} for the simulated data set to a value of 200 keV. In this case, the constraints on E_cut are recovered with regular monomodal posteriors.

  Γ–E_cut dependence: For models which have E_cut, Γ shows a dependency on E_cut when it is kept free, this is clear from the skewness in Γ-distribution when E_cut is a variable parameter in the fits. For BORUS when SCPL-ICPL are uncoupled, posteriors of Γ in the case of E_cut frozen fits are more asymmetric and have higher values of R₉₀ and δ. For UXCLUMPY, the same situation holds, but the availability and dependence on SCPL down to 0.3 keV in XMM–Newton, diminishes the effect of the E_cut in the fit where SCPL and ICPL are coupled. The overall trends suggest that E_cut might influence constraints on ICPL-Γ.

• Angle of Inclination (θ_i): The angle of inclination (θ_i) to the torus axis is determined the flux and the shape of the scattered continuum reaching the observer. Thus, the models will use the CRH and/or the low energy (E < 6.0 keV) tail of the scattered continuum, to estimate θ_i. For UXCLUMPY, in the UX-OMNI setup, in the XMM–Newton+NuSTAR joint fits, the value of δ for θ_i is 0.41, indicating that it is not a well constrained parameter. For CTORUS, θ_i is skewed towards the higher angles. This is an artefact of the correlation of θ_i with N_cloud (Fig. 6a). The N_cloud posterior which covers the whole allowed range (see pt. vii) causes the θ_i posteriors to get skewed to higher angles. In models based on bi-conical cutouts like BORUS and ETORUS θ_i is constrained well if the torus opening angle (⁠|$\rm cos \theta _{\rm o}$| or equivalently C_F,tor of BORUS) constrained well. For cases when θ_o shows trends of wider posteriors or bimodality θ_i shows the same trend. Additionally bimodality in θ_i was seen if θ_o shows bimodality. In doughnut-based models like MYTORUS and RXTORUS, N_{H, eq} is the input parameter and N_{H, los} = f(N_{H, eq}, θ_i, c/a). Because of the dependence N_{H, los} is a well constrained parameter in this regime and thus can help narrow down the favourable θ_i region in the parameter space. Further localization and limiting the N_{H, eq} – c/a – θ_i degeneracy is provided by features of the scattered continuum. We find that the N_{H, eq} – θ_i correlation is strong when c/a is frozen. However, when c/a is kept free (in RXTORUS), the N_{H, eq}–θ_i correlation weakens and the c/a – θ_i correlation strengthens.

• Relative normalization between the transmitted and the reflected components (T/R): In many models, the normalization of the zeroth-order continuum and the scattered continuum can was kept equal i.e. T/R = 1 or kept free as discussed in Section 2.2. We simulated data for two cases, (A) |$T/R\, \gt 1$| for CTORUS (T/R = 1.8) and MYTORUS (T/R = 1.5) (B) T/R = 1 for all models. We kept T/R free for only MYTORUS, CTORUS, and ETORUS. For all these models, in both case-A and case-B the parameters were recovered. For MYTORUS, CTORUS, and ETORUS (r, δ) were found to be [A:(0.77, 0.08), B:(0.12, 0.06)], [A:(8 × 10⁻³, 0.154), B:(0.09, 0.17)], [B:(0.9, 0.29)], respectively.

• Relative Normalization of the SCPL (C_sc-pl): For the models MYTORUS, RXTORUS, BORUS, ETORUS, and CTORUS C_sc-pl are recovered with δ ∼ 0.067–0.54. For UXCLUMPY, in the UX-OMNI, C_sc-pl posteriors were recovered with values of δ and r equal to 0.54 and 0.28, respectively. C_sc-pl is strongly correlated with σ_o, θ_i, and C_frac of the inner ring. In the UX-ZCPL setup the torus C_sc-pl is recovered with δ ≃ 0.1, C_sc-pl is correlated with σ_o and θ_i, but not with C_frac.

• Parameters of Torus morphology:

  • c/a of RXTORUS: The posteriors returned from the fits are regular monomodal distributions. The c/a-parameter shows strong correlation with the posteriors of N_{H, eq} and θ_i. c/a in RXTORUS determines the opening angle of the doughnut, hence the c/a – θ_i correlation. The constraints on c/a are tight as δ ≃ 0.086 and r ≃ 0.024.

  • Opening angle (θ_o): For ETORUS θ_o is a direct input. However, in BORUS the opening angle is given by θ_o = cos ⁻¹C_{F, tor} for a contiguous torus. θ_o shows different trends in ETORUS and BORUS. For ETORUS, θ_o returned very broad posteriors with δ = 0.2 when T/R was frozen to the input value (T/R = 1) and worsened to 0.35 when T/R was free, with strong θ_o–T/R correlation. For BORUS (T/R = 1 frozen in all cases) however the constraints are better with δ ≃ 0.02. The parameter constraints are determined from the shape of the lower energy (E < 6.0 keV) tail of the reflected continuum and the shape of the CRH.

  • Tor-sigma (σ_o) of UXCLUMPY: When data simulated in the UX-OMNI setup and was fit it with the same, σ_o is recovered with δ ≃ 0.33. The 90 per cent confidence region was spread over almost |$\sim 65{{\ \rm per\ cent}}$| of the prior range for the given flux level (see Fig. 3). The most distinct influence σ_o and has is the soft band tail of the reflection component. σ_o along with C_frac also affects the relative height of the CRH to the FeKα absorption edge. Additionally, the zeroth-order continuum is strongly present in the 7–25 keV band and dominates the reflection component. This invokes a strong σ_o–C_frac–C_scpl correlation and a large spread in the posteriors of σ_o (see Fig. 4 or 3). When the UX-OMNI data set is fit with the UX-ZCPL setup, σ_o was not recovered. While we get a lower value of δ ≃ 0.2, the value of r increases to 2.2, indicating a discrepancy. The SCPL excess was adjusted by decreasing σ_o, which increased the amount of the scattered component to make up for the comparatively higher contribution from the warm mirror. Fig. 4 shows the comparison of the UXCLUMPY morphological parameters, for the UX-OMNI and UX-ZCPL cases.

  • Covering fraction (C_frac) of inner ring of UXCLUMPY: In the UX-OMNI setup, there exists a strong C_frac–σ_o degeneracy (see Fig. 3). The posteriors were recovered with δ ≃ 0.14. But strong degeneracy with σ_o can result in discrepancy with the input. The constraints on the C_frac are derived from the energy band of the scattered continuum between 7 and 25 keV. In the MCT regime, the dominance of the zeroth-order continuum over the scattered continuum in the 7–30 keV energy band results in poor constraints on C_frac. In the UX-ZCPL setup, C_frac was recovered with δ ≃ 0.08. From the contour plots (see Fig. 4), it can be seen that the nature of σ_o–C_frac correlation is different for the two setups. UX-OMNI shows a σ_o–C_frac Pearson’s correlation coefficient of 0.75 (N_samples ≃ 3718), whereas for UX-ZCPL setup we get σ_o–C_frac correlation of ∼0.14 (N_samples ≃ 2682) with a confidence of greater than 99.99 per cent.

  • N_cloud of CTORUS: The input value of N_cloud of data simulated with CTORUS was set to 4. We simulated the data sets for |$\rm \mathit{ T}/\mathit{ R} = 1.0$| and 1.8, respectively (see point no. iv). The values of R₉₀ and δ for both the spectra with |$\rm \mathit{ T}/\mathit{ R}=1.8$| and T/R = 1.0 are similar (R₉₀ ∼ 5 and δ ∼ 0.6). The constraints on N_cloud are determined from two features of the reflected continuum: the E < 6 keV tail and the CRH. The reflected continuum has lower flux by a factor of ∼1/5 (T/R = 1.8) to ∼1/3 (T/R = 1) than the zeroth-order continuum in the 3–100 keV band in this regime; constraints from CRH become overwhelmed, leading to flattened posteriors of N_cloud. Additionally, the soft tail (E < 6 keV) of the scattered component is weaker compared to other models, leading the SCPL to be more dominant in this band, weakening the constraints on N_cloud.

  • Iron abundance (A_Fe) of BORUS: The zeroth-order continuum is simulated using the cutoff power law absorbed through zTBABS×CABS and abundance is fixed at 1. Thus, technically there is no way to consistently simulate a generic spectrum with any value of A_Fe by linking the zeroth-order continuum-A_Fe to the scattered component-A_Fe. Thus, to make the spectrum consistent, for the data simulated under BORUS: A_Fe = 1.0 for the scattered component. A_Fe was recovered (r_q ≃ 0.2) with δ ≃ 0.05 by measuring the features e.g. depth of the Fe K absorption edge, the height of the FeKα line of only the reflected component.

• apec components: apec parameters were genereally recovered when the lowest energy bin of the data is lower than the apec temperature included in the spectrum. Thus with this knowledge we kept the apec temperatures and normalizations frozen at the input values, except for a medium Compton-thick fit of UXCLUMPY where all the apec parameters are kept free.

Figure 3.

Contours for UXCLUMPY IM analysis in the MCT regime, with N_{H, los} = 100 as input. The degeneracy or correlation between the parameters (logC_sc-pl–σ_o–C_frac) with extremely skewed contours can be noticed. The lines or cross-hairs coloured red denote the input value and those coloured green denote the median value calculated from the posterior distribution.

Figure 4.

Contour plots indicating the correlation of the parameters C_sc-pl (relative normalization of warm-mirror or scattered power law), σ_o, C_frac, and θ_i of UXCLUMPY, when data simulated under UX-OMNI is fit with the (a) UX-OMNI (b) UX-ZCPL setup. The red and the green lines in the histograms and cross-hairs in the contour plots mark the input and the median values calculated from the distributions, respectively. The nature of the correlation between the parameters changes depending on the SCPL model. In this case we keep the temperatures and normalizations of the apec components and the instrumental constants frozen at their input values.

3.2 Heavy Compton-thick regime

The photon-index Γ did not show significant changes in the constraints. However several other parameters showed a significant change in terms of recovery, degeneracy, and constraints in the HCT regime compared to the MCT regime. We summarize the list of values of δ and r for all parameters and models for the MCT and HCT regime in Table 2. We discuss the results of fits for data simulated in the HCT regime, with N_{H, los} ≃ 500 and 2–10 keV flux values similar to that of the MCT for all models. The most significant difference of this regime is the zeroth-order continuum gets overwhelmed by the scattered continuum (⁠|$\left[ \frac{F_{\rm \mathit{ T}}}{F_{\rm \mathit{ R}}}\right]_{\rm 3-100\, keV} \simeq 0.1$|⁠). We discuss the results from the HCT regime in this subsection.

• Column density (N_H): The constraints on N_{H, los} worsen (top panel of Figs 7 and 8a) in this regime because of the decrement of |$\left[ \frac{F_{\rm \mathit{ T}}}{F_{\rm \mathit{ R}}}\right]_{\rm 3-100 \,keV}$| by an order of magnitude. δ ranges between 0.03 and 0.21 (average δ ≃ 0.09, table 2) in the HCT-regime for IM fits of all models. For estimating the constraints on N_{H, los} in the HCT regime (in the absence of the zeroth-order continuum) the features of the scattered continuum e.g. the relative height of the CRH to the Fe K-edge, the E < 5 keV tail and the shape of the Compton hump play the most significant role.

• Angle of inclination (θ_i): In all cases (except ETORUS), the constraints in θ_i either improved or remained same as the MCT in the HCT regime. This improvement can be attributed to the improved prominence of the scattered continuum resulting in better constraints in the morphological parameters which affects θ_i posteriors. For example in CTORUS, posteriors of N_cloud in the HCT1 case improved, thus improving δ for θ_i to 0.05 (Fig. 6). For UXCLUMPY the parameter posteriors of σ_○ and C_frac improved, leading to an improvement in the value of δ for θ_i. Although for the case of RXTORUS and MYTORUS, the constraints do not improve but remained almost the same as that of the MCT regime. However, ETORUS showed worse constraints on θ_i, when T/R is free but better when T/R if frozen (see next bullet point on θ_o). In the HCT regime, as an average generic trend posterior on θ_i are better constrained, with an average value of δ lower than the MCT regime for the respective models (Table 2).

• Parameters of morphology:

  • c/a of RXTORUS: The constraints on c/a improved (Fig. 8b) in this regime. The parameter was recovered, with δ ≃ 0.04 which is nearly half of that calculated in the MCT regime. The value of r_q on the other hand became worse (increased to 0.63 from 0.02 in the MCT regime), (see Table 2) compared to the MCT regime. One reason for this change is the decrease in the spread of 90 per cent confidence region (R₉₀) in the HCT regime, thus increasing the ratio (Δq/R₉₀).

  • Opening angle (θ_o) or cos θ_o: For ETORUS the values of δ improved mildly to 0.14 (from 0.2 in the MCT) when T/R is kept frozen at input in the HCT regime. When T/R is kept free, the constraints worsen (δ ≃ 0.8) compared to MCT and are spread over the whole prior range. The soft tail (E < 6 keV) of the scattered continuum of ETORUS is considerably lower compared to the other models in the HCT regime. For BORUS we observe bimodality, with two disjoint peaks in the posteriors. The major peak corresponds to the input and the minor peak is significantly discrepant from the input. We can evaluate the quantiles independently for each peaks. When compared to the MCT regime the values of R₉₀ and δ are similar. R₉₀=1.51 and 2.12 for the major and minor peak, respectively. The ratio of the number of samples corresponding to major and minor is N_M/N_m ≃ 17. It may be possible to reject the minor peak as it corresponds to lower probalility. The bimodality in θ_o induces bimodality in θ_i and C_sc-pl.

  • Tor-sigma (σ_o) of UXCLUMPY: Because improved prominance of the scattered-continuum in the HCT regime, the influence of the effects of σ_○ on 6–40 keV becomes evident unlike in MCT. In the HCT when the σ_o decreases, it increases the level of the soft band (E < 6 keV) tail, it increase the warm-mirror component if the spectra are simulated in the UX-OMNI setup and it approaches the Compton-thick ring and that introduces a mild depression in the 7 and 30 keV band similar to that of the Compton-thick inner ring. When σ_○ approaches 90^○ the system approaches an almost 4π covered system, diminishing the soft band significantly. These features can potentially establish stronger constraints as observed (comparison of Fig. 5 with 3 and posterior overplot in Fig. 8d) in the HCT compared to MCT as the scattered continuum is the dominant part of the spectrum.

  • Covering fraction (C_frac) of inner ring of UXCLUMPY: The constraints improve (Fig. 8c) because of improved prominance of the features of the scattered continuum. This improved the posteriors of σ_○ and C_sc-pl, leading to elimination of the heavy correlation seen otherwise in the MCT regime. We found that δ for C_frac improves from 0.14 in the MCT regime to 0.08 in the HCT regime.

  • N_cloud of CTORUS: For the HCT regime we test two cases:

    1. The SCPL is stronger and contributes much more flux (HCT0) than the scattered continuum from the torus (green curve in Fig. 2d) in the 2–50 keV band [f_torus/f_total]₂₋₁₀ = 0.45. This case represents objects with higher diffused/scattered emission. For this case, we get δ = 0.5 and R₉₀ ≃ 4.18.

    2. The SCPL is weaker than the soft tail of the scattered continuum from the torus (HCT1, magenta curve in Fig. 2d): [f_torus/f_total]₂₋₁₀ ≃ 0.90. In this case, the CRH is much stronger and represents cases where the diffused or scattered emission is weak but the torus continuum is strong. For this case we get δ = 0.2 and R₉₀ ≃ 1.97. This is because of the presence of the strong soft-tail from the torus scattered continuum.

    For both cases, the posteriors of N_cloud (Figs 6b and 8e) improve in the HCT regime, as the CRH and E < 4 keV tail become dominant over the zeroth-order continuum.

  • Iron abundance (A_Fe) of BORUS: The constraints of BORUS improved in the HCT regime (Fig. 8f), δ ≃ 0.02. The problems introduced because of the inconsistency of the iron abundance in the zeroth-order continuum and the scattered continuum would be minimal or absent in the HCT regime, because of the low value of |$\left[ \frac{F_{\rm \mathit{ T}}}{F_{\rm \mathit{ R}}}\right]_{\rm 3-100\, keV}$|⁠, which almost eliminates the zeroth-order continuum, thus enabling us to calculate A_Fe from the features of the dominant scattered continuum only. We found that δ improves to 0.02 in the HCT-regime from 0.07 in the MCT regime (see Table 2).

Figure 5.

Contour plot for UXCLUMPY IM analysis in the HCT regime, with N_{H, los} = 500 as the input. A noticeable difference exists in the parameter posteriors in terms of (1) better posterior localization/constraints on parameters of morphology and (2) lesser level of parameter degeneracy exhibited by the nature of the 2D contours compared to the MCT regime (Fig. 3). The lines or cross-hairs coloured red denote the input value and those coloured green denote the median value calculated from the posterior distribution.

Figure 6.

Contours for CTORUS IM analysis in the (a) MCT: input N_{H, los} = 100, (b) HCT0: input N_H = 500, (b) HCT1: input N_{H, los} = 500. The lines or cross-hairs coloured red denote the input value and those coloured green denote the median value calculated from the posterior distribution. The influence on N_cloud posterior on θ_i due to its poor localization is noticeable in (a). But constraints improve considerably in (b).

Figure 7.

We plot the posterior medians and 90 per cent confidence ranges for the parameters that are common across all the models. The diamond markers are for the HCT regime and the square markers for MCT. Noticeable difference in the trends of parameters across regimes for different models can be seen e.g. N_{H, los} is well constrained in the MCT regime for all models. However, the extent of worsening of the constraints on N_{H, los} in the HCT regime is different for different models.

Figure 8.

Overplot of the distributions of X_norm = X_posterior/X_input for MCT and HCT regime, where X_posterior represents the point corresponding to the posterior distribution and X_input is the input value used for data simulation. Noticeable differences exist in the distribution of the selected parameters in the MCT and HCT regimes. e.g. subplot (a) shows that the probability distribution of N_H(note: in linear space) is more localized in the MCT regime. Posterior distribution of other parameters get more localized in the HCT regime by virtue of higher flux contribution from the scattered continuum, as shown in subplots (b)–(f).

Figure 9.

We compare the 1D posterior distributions calculated from XMM–Newton, NuSTAR, and joint XMM–Newton + NuSTAR fits in this figure. We notice significant dependence in the nature of the posterior distributions on the instrument and their combinations. These posteriors prove that for decent constraints on some parameters e.g. C_frac of the inner ring of UXCLUMPY, NuSTAR data is a requirement, whereas for N_cloud in the HCT regime joint data from XMM–Newton and NuSTAR is required.

Figure 10.

Contour plot results for the analysis M_CTOR → D_MYT in MCT regime. The posterior of N_cloud is irregular as there exists no maxima in the allowed prior range, but a minima, with ridges at the prior edges. This ‘bifurcation’ in the PDF of N_cloud gives rise to bimodality in the other parameters.

Figure 11.

(a) Contour plot results for the analysis M_RXT → D_MYT in MCT regime. The noticeable discrepancy can be seen in the parameter c/a which in turn influences the constraints on logC_SCPL, θ_i posteriors. (b) spectral overplot of the input model (green), the best-fitting spectrum from M_RXT → D_MYT (blue) and the RXTORUS model under MYTORUS input parameters (red). The dashed spectra are the scattered continua of the respective model spectra. The scattered components of RXTORUS and MYTORUS spectra under same parameters are different because of the difference in radiative process assumptions. χ²/dof = 1.11 for the best-fitting parameters.

Figure 12.

We show the behaviour of morphological parameters of UXCLUMPY for M_UXCL → D_j-fits in the MCT regime(a)σ_o and (b)C_frac of the inner ring. σ_o is always ≤30^○ and C_frac in most cases have irregular posteriors consistent with 0 and a regular posterior with median value consistent ≃ 0.1. The blue markers with errorbars represent the median values and uncertainty on σ_o and C_frac for typical of exposure (diamond marker: t_XMM = 100 ks, t_NuSTAR = 50 ks) and a deeper observation (circle marker: t_XMM = 1.0 Ms, t_NuStar = 0.5 Ms) for M_UXCL → D_MYT. The dashed vertical lines mark the median values obtained from the spectrum with high exposure. The median values and the errors obtained from the typical and deeper observations show that trends are systematic and occurring due to model difference. See discussion in Section 4.1.

Figure 13.

Cross- and intra-model fits for the case M_j → D_UXCL in the MCT regime, j = MYT, RXT, CTOR, and BOR models. The discrepancies in Γ are typically rather large, with a flattening by 0.2–0.3, and potentially impacting conclusions about the coronal power law. In contrast, discrepancies in N_{H, los} are much smaller, at only |$\sim 10 {{\ \rm per\ cent}}$| which does not significantly impact the spectral fit.

Table 2.

In this table, we compare the average values of δ and r_q for selected parameters in the MCT and HCT regimes. The averages are taken over all analysed sub cases of all intramodel fits of the models that have the given parameter. For most parameters of morphology, we get a better constraints except θ_o of ETORUS which worsens and cosθ_o of BORUS which does not show any significant improvement in the HCT regime.

Parameter	δ		rq
	MCT	HCT	MCT	HCT
NH, los	0.02	0.09	0.44	0.51
Γ	0.018	0.014	0.35	0.46
θi	0.18	0.13	0.60	0.49
cosθo of BORUS	0.035	0.021	0.84	0.23
θo of ETORUS	0.31	0.47	0.5	0.67
σo of UXCLUMPY	0.32	0.2	0.37	0.37
Cfrac of UXCLUMPY ring	0.13	0.08	0.51	0.84
Ncloud of CTORUS	0.64	0.26	0.47	0.8
AFe of BORUS	0.059	0.026	0.41	0.49
c/a of RXTORUS	0.086	0.043	0.024	0.63

3.3 Bulk simulations

The properties of the posterior distributions e.g. R₉₀, δ, and r_q are affected by statistical processes involved in data generation. This implies that for a given theoretical spectral model, the data sets generated in each case will be different and Bayesian analysis performed on these data sets will return posteriors of parameters with median values, R₉₀, δ, and r_q different from each other. However these distributions will follow a characteristic distribution for each parameter. To verify this we perform bulk simulations of spectra for the UXCLUMPY in the MCT and HCT regime. We limit ourselves to simulating 100 spectrum for each case, due to the limits of computational power and time. We perform intramodel fits on each of these spectra and calculate posteriors in each case.

We calculate the distribution of certain properties of the posteriors (e.g. R₉₀) of different regimes. We find several properties of the distribution which are consistent with our observations from the single fits. N_{H, los} remains well constrained in the MCT regime compared to the HCT regime. The distributions of R₉₀ and δ determined from the bulk fitting of the spectra of the torus-morphological parameters are localized to the lower values for the HCT regime compared to those of the MCT regime which is spread towards comparatively higher values. However, the distribution of R₉₀ from the two regimes shows overlap with each other as shown in Table 3. Thus the values of R₉₀ and δ estimated in the single spectral fit in the previous sections are parts of these distributions of respective parameters. That is for a given parameter the probability that the value of R₉₀ (or δ) calculated from a single fit will remain enclosed in the respective values and the errorbar on the parameter quoted in Table 3 is 90 per cent. The implies that δ values quoted in Table 2 which reflect the statistical errors, should not be considered as absolute, as they can be affected by the statistical process involved in the data generation, implying a chance of getting slightly different values of R₉₀ and δ for data sets simulated for the same model under same conditions. The extent of variation in the R₉₀ and δ depends on the regime and the parameter of interest. Additionally, Table 4 shows that a decrease in R₉₀ or δ values results in a decrease in the uncertainty on the same. In the context of intramodel fits, this means, as a parameter becomes better constrained, the statistical error on that parameter becomes less influenced by statistical processes involved in data simulation or generation.

Table 3.

In this table, we show the distribution of the parameters characterizing posteriors distributions, calculated from the bulk fitting of 100 spectrum. The ‘overlap percentage’ (abbreviated as per cent-overlap) is the fraction of the values of R₉₀ or δ from the wider distribution that overlaps with values lesser than the 0.95th quantile of the narrower distribution. The narrower which for the Γ and N_{H, los} is MCT regime and for σ_o, C_frac, and θ_i is the HCT regime. These values give an idea about how much the statistical errors (R₉₀, δ, and r_q) might vary due to random process associated with the data generation.

Parameter	R90			δ			rq
	MCT	HCT	Per cent-overlap	MCT	HCT	Per cent-overlap	MCT	HCT
Γ	|$0.080^{+0.013}_{-0.015}$|	|$0.081^{+0.018}_{-0.017}$|	∼100 per cent	|$0.021^{+0.004}_{-0.004}$|	|$0.022^{+0.005}_{-0.005}$|	∼100 per cent	|$0.40_{-0.38}^{+0.85}$|	|$0.63_{-0.50}^{+0.86}$|
NH, los	|$3.24^{+0.29}_{-0.31}$|	|$143^{+251}_{-57}$|	∼0 per cent	|$0.016^{+0.0016}_{-0.0017}$|	|$0.138^{+0.211}_{-0.048}$|	∼0 per cent	|$0.56_{-0.52}^{+0.69}$|	|$0.32_{-0.28}^{+0.76}$|
σo	|$26^{+8}_{-6}$|	|$20^{+5}_{-4}$|	43 per cent	|$0.27^{+0.09}_{-0.08}$|	|$0.22^{+0.03}_{-0.03}$|	38 per cent	|$0.38_{-0.34}^{+0.91}$|	|$0.49_{-0.43}^{+0.82}$|
Cfrac
of UXCLUMPY ring	|$0.094^{+0.039}_{-0.016}$|	|$0.07^{+0.008}_{-0.013}$|	10 per cent	|$0.11^{+0.04}_{-0.02}$|	|$0.09^{+0.01}_{-0.02}$|	32 per cent	|$0.53_{-0.48}^{+0.66}$|	|$0.47^{-0.44}_{+0.73}$|
θi	|$52^{+13}_{-14}$|	|$29^{+9}_{-7}$|	6 per cent	|$0.47^{+0.60}_{-0.19}$|	|$0.28^{+0.2}_{-0.1}$|	53 per cent	|$0.41_{-0.40}^{+0.78}$|	|$0.53_{-0.45}^{+0.90}$|

Table 4.

We summarize the input parameters that are used for data simulation under UXCLUMPY and the median values and errors calculated from the parameter posteriors in this table, for the cases M_j → D_UXCL, where j = UXCL, MYT, RXT, CTOR, and BOR in the MCT regime. For parameter estimation the full energy range allowed by a model has been taken into consideration and the χ²/dof value is also calculated over the entire energy range. The values of evidence and Bayes-Factor are quoted for the case when all the energy range of the data used for analysis is 1.2–78.0 keV for all models.

	Input	Fitting model
Parameters	Input for DUXCL	MUXCL	MMYT	MCTOR	MRXT	MBOR (coupled)
Csc-pl	−2	|$-1.93^{+0.24}_{-0.18}$|	|$-2.28^{+0.04}_{-0.04}$|	|$-2.19^{+0.05}_{-0.09}$|	|$-2.24^{+0.05}_{-0.05}$|	|$-2.49^{+0.03}_{-0.04}$|
NH, eq	–	–	|$211^{+40}_{-34}$|	–	|$218^{+20}_{-25}$|	–
NH, los	100	|$100^{+1.7}_{-1.7}$|	92	|$104.5^{+0.24}_{-0.29}$|	92	|$91.83^{+1.5}_{-1.6}$|
Γ	1.9	|$1.88^{+0.04}_{-0.05}$|	|$1.72^{+0.03}_{-0.03}$|	|$1.71^{+0.03}_{-0.03}$|	|$1.8^{+0.03}_{-0.03}$|	|$1.60^{+0.02}_{-0.02}$|
Ecut	400	|$323^{+67}_{92}$|	–	–	–	500 (frozen)
σo	45	|$49^{+16}_{-17}$|	–	–	–	–
Cfrac	0.4	|$0.42^{+0.07}_{-0.05}$|	–	–	–	–
Ncloud	–	–	–	|$3.36^{+3.95}_{-1.2}$|	–	–
T/R	1 (frozen)	1 (frozen)	|$1.90^{+0.10}_{-0.10}$|	|$3.03^{+0.53}_{-0.44}$|	|$2.42^{+0.16}_{-0.15}$|	|$6.38^{+0.73}_{-0.73}$|
c/a	–	–	–	–	|$0.88^{+0.03}_{+0.02}$|	–
CFrac, tor or |$\rm cos \theta _{\rm o}$|	–	–	–	–	–	|$0.85^{+0.06}_{-0.12}$|
θi or |$\rm cos \theta _{\rm i}$| (for BOR and CTOR)	60.0	|$64.56^{+20}_{-33}$|	|$63.25^{+1.43}_{-0.92}$|	|$0.53^{+0.29}_{-0.04}$|	|$36.85^{+3.48}_{-2.62}$|	|$0.88^{+0.05}_{-0.14}$|
AFe (for BOR)	–	–	–	–	–	|$0.96^{+0.06}_{-0.09}$|
L7–50/Lin, 7–50 of coronal power law	1.0	1.04 ± 0.04	0.87 ± 0.02	1.08 ± 0.04	0.89 ± 0.02	0.96 ± 0.03
χ2/dof	–	1.004	1.14	1.09	1.11	1.34
For evidence estimates: 1.2 ≤ E ≤ 78 keV
|$\rm log\mathit{ Z}$|	–	−592.96	−626.01	−654.60	−625.09	−800.12
|$\rm log BF = log \mathit{ Z}_{j}-log \mathit{ Z}_{UXCL}$|	–	0	−33.06	−61.64	−32.13	−207.16

Finally, as a test of validation of the errors for individual fits, we compare the average value of R₉₀ (〈R₉₀〉 = Σ_i(q_0.95 − q_0.05)) over all the individual fits (given a parameter) and compare it with the extent of 90 per cent confidence region of the distribution of median values of the posteriors (R_{90, med}). We find that values of 〈R₉₀〉 and values of R_{90, med} are very similar for both regimes. For example, values of 〈R₉₀〉 and R_{90, med} for Γ are ∼0.81 and 0.60 for MCT regime and 0.81 and 0.64 for HCT regime, respectively. The values of 〈R₉₀〉 and R_{90, med} for σ_o are, respectively, 26.14 and 20.55 for the MCT regime, and 20.04 and 18.68 for the HCT regime. These similarities indicate the properties fo the uncertainties obtained from individual fits are valid in that they are representative of the distributions of medians.

3.4 Instrument dependence of model parameters

In this section, we seek to understand which parameter constraints are leveraged by which instrument, and hence understand the required energy band coverage for measurement of a particular parameter. We fitted either XMM–Newton data only, or NuSTAR data only, and compare the results with those obtained above from joint fits. This was done for both MCT and HCT scenarios. We summarize our primary findings as follows:

• Column density (N_{H, los}): In the MCT regime N_{H, los} parameter is more reliably determined with the XMM–Newton data which facilitates the better constraints compared to NuSTAR because of better energy resolution, thus facilitating better measurement of the Fe K absorption edge and the rollover at E < 5 keV. In the HCT regime the scattered component is dominant and the strong attenuation of the transmitted component does not allow the measurement of a rollover or edge of the transmitted component. Both the instruments individually do not provide reliable measurement in the HCT regime. However joint fits in all cases are more reliable in terms of values returned and in terms of low statistical errors, which makes joint observation preferable when proper measurement of N_{H, los} is concerned.

• σ_o of UXCLUMPY: σ_o of the cloud distribution is not a well constrained parameter in any case, with values of δ as large as ∼0.30. The posteriors are wider (Fig. 9d) for NuSTAR compared to XMM–Newton only, as the E < 4 keV tail of the scattered component is absent in NuSTAR only data.

• Covering fraction (C_frac) of inner ring of UXCLUMPY: The constraints are generally drawn from the shape of the CRH, that is, 7 < E ≲ 30 keV band. We found that fitting only XMM–Newton data returns very wide irregular posterior, whereas NuSTAR only returns regular monomodal distribution, with stronger constraints. The constraints improved only by a small amount in the joint analysis. Thus, data from NuSTAR are a compulsory requirement to constrain the covering fraction (C_frac) of such a ring if present.

• c/a of RXTORUS: The values of δ are similar for both the instruments. But the median value of NuSTAR fit is discrepant compared to the input and posteriors are wider whereas the XMM–Newton fit is not. On the other hand, the constraints on the lower limit, are stronger in the XMM–Newton fit, whereas constraints on the higher limit are stronger for NuSTAR fit. The constraints on the joint fits, thus depend strongly on data from both the XMM–Newton and NuSTAR instruments. The median value of NuSTAR posteriors differ from the input by nearly 24 per cent, due to unavailability of the E < 4 keV tail of the scattered component.

• N_cloud of CTORUS: XMM–Newton-only fits provides tighter constraints compared to NuSTAR as shown in the figure (e.g.R_{90, Nu} = 6.48 and R_{90, XMM} = 4.35 for HCT). The joint posteriors improve significantly (e.g. R_{90, joint} = 1.97 for HCT) under the influence of XMM–Newton. Thus, for a more localized posterior on N_cloud, joint fits with XMM–Newton data is a requirement.

• A_Fe of BORUS: Both the XMM–Newton and NuSTAR parameters have values quite different from the joint fits and the input. The posteriors of NuSTAR are wider. We would require both NuSTAR and XMM data to constrain the posteriors better accuracy.

Our analysis of the data from the instruments separately, also shows that simultaneous coverage of 0.3 < E < 78 keV band improves the posteriors significantly, compared to separate E ∼ 10 keV through XMM–Newton or E > 4 keV only through NuSTAR coverages. We jointly show the 1D posterior plots of XMM–Newton only, NuSTAR only and XMM–Newton plus NuSTAR joint fits for a few parameters in Fig. 9.

4 RESULTS: CROSS-MODEL FITS

In this section, we discuss the results of cross-model fits (CM-fits hereafter). We use the notation M_j → D_i to denote that the model j is fit to data simulated under model i. The indices i and j are the abbreviations for the model of which the reflection component is featured in the data simulations or fits. The abbreviations for the torus models are summarized in Table 1. For example, when the fitting model is UXCLUMPY, the notation would be M_UXCL → D_MYT. For other cases, the abbreviations for model combinations are defined locally. The CM fits will identify and estimate the level of degeneracy between the models, show which parameters adjust themselves to fit a given data set when the assumed morphology is wrong. We do not test all possible model combinations; instead, we concentrate on only the following cases:

• M_MYT → D_j, |$\rm j =$| RXT, UXCL: spectral change due to difference only in radiative process and how well a simple smooth doughnut can represent complex arrangements.

• M_RXT → D_j, |$\rm j =$| MYT, UXCL: similar to MYT but for the variable doughnut parameter (c/a).

• M_BOR → D_j, |$\rm j =$| CTOR, UXCL: what happens when a model with a simplified contiguous cut-out geometry is applied to clumpy morphologies.

• M_CTOR → D_j, |$\rm j =$| MYT, UXCL, BOR: how a clumpy torus model with an estimate on the number of clouds in the radial direction acts when applied to data simulated under other contiguous and clumpy models.

• M_UXCL → D_j, |$\rm j =$| MYT, RXT, BOR, CTOR: same as CTOR but UXCLUMPY has a more complex cloud distribution.

Different models treat and make predictions for the shape/intensities of the emission complex of Fe Kα and other fluorescent emission lines due to different scattering geometries and different physical processes in their radiative transfer codes e.g. MYTORUS contains Fe-lines with Compton shoulders only and other models contain a wider range of species. To account for such discrepancy between models, one can perform cross-model fitting ignoring certain narrow energy ranges thus focusing more on constraints from the continuum (see Section 4.2). We can expect that there may be narrow ranges of data-model residuals corresponding to differently modelled Fe Kα and other line complexes. Future work on model comparison can be done when models offer a more consistent treatment of narrow features like fluorescent lines and Compton shoulders.

4.1 Medium Compton-thick regime

• Column density (N_{H, los}): Across all models tested in the MCT regime, LOS column density (N_{H, los}) is the only parameter that was practically recovered. The dominance of the zeroth-order continuum in this regime is the reason for the well-determined constraint and values consistent with the input for N_{H, los} irrespective of the model. The maximum discrepancy was observed for M_CTOR → M_BOR, where the percentage discrepancy corresponds to 31.8 per cent and r_q corresponds to 8.2. With regard to specific trends we find that the clumpy torus models returned higher values of N_{H, los} when fit to contiguous torus and the converse also holds true. The values of r_q correspond to the variation of N_{H, los} by 30 per cent. In order to check the discrepancy statistically we also perform bulk CM-fits of the case M_MYT → D_UXCL to 100 spectra simulated under UXCLUMPY. The minimum and maximum discrepancy in N_{H, los} in the bulk fits were found to be 1.7 and 8.6 per cent, respectively. This discrepancy in the values of N_{H, los} will not create profound physical inconsistencies and will not significantly change the scientific implications in terms of the spectral shape and the components present. Thus, we conclude that N_{H, los} returned from the fits were mostly consistent with the input in the MCT regime.

• Photon-index (Γ): Both the CRH and zeroth-order continuum determine the value of Γ. The values of r_q that we observed for Γ in the CM-fits range from 0.67 to 15 with Γ_in = 1.90 for all tested cases. The best value was for the case M_MYT → D_RXT (Γ_out = 1.91) and the worst value was obtained for M_BOR → D_UXCL. In the MCT regime, for a specific case M_j → D_UXCL, where j = MYT, RXT, BOR, and CTOR, the fits systematically returned flatter values of Γ, 1.60 ≲ Γ_out ≲ 1.75 (see Table 4). In the converse case viz. M_UXCL → D_j, where j = MYT, RXT, BOR, and CTOR, the estimated values of Γ were in the range 2.02 ≲ Γ_out ≲ 2.17. The high discrepancy (quantified by r_q) in values of Γ is a direct implication of low statistical errors, caused by the intrinsic difference in the CRH shape for different models. Therefore, we can conclude that the retrieved value of Γ is thus strongly model dependent (see Fig. 13).

• Angle of inclination (θ_i): The values of θ_i returned from the CM fits are dependent on the torus morphology. However, there might be cases where interpretation of the values can be contextual and is dependent on morphology. We observed specific trends for certain models in the CM-fits. We list these in the following bullet points:

  • In CTORUS the zeroth-order continuum is independent of θ_i. The information about the only parameters of morphology (N_cloud) and geometry (θ_i) is in the scattered component. In CTORUS, the obscuring clouds are located in a hypothetical biconically cutout thick shell between 60^○ and 90^○ (Liu & Li 2014) with respect to the symmetry axis. In the M_CTOR → D_UXCL fit when we assume a prior [0^○, 90^○], we get an inconsistent value of θ_i < 60^○ according to the geometry proposed in Liu & Li (2014). However, a possible misinterpretation is the existence of a stray clump(s) or obscurer at θ = θ_i < 60^○ which is outside the hypothetical shell enclosure of CTORUS. This stray clump or obscurer can be thought to be contributing to N_{H, los}. This is just a result of spectral or model difference and hence an artefact of fitting (also see the similar discussion for cos θ_o for BORUS). When we limit the prior range to [60^○, 90^○], the posterior monotonically increases towards the 60^○ edge.

  • For the CM-fits, M_UXCL → D_j where j = CTOR, MYT, BOR, θ_i assumes an edge-on configuration with monotonically increasing posteriors consistent with 90^○. The only exception is M_UXCL → D_RXT, where the median value calculated from the posterior is ≃ 66^○ for an input θ_i = 70^○; we suspect that the excess contribution in the softer bands from Rayleigh scattering might be the result of the lower values of θ_i. The input value of θ_i and the value returned from fit M_UXCL → D_RXT value might agree by chance.

  We discuss the case of θ_i for the M_RXT → D_MYT case along with c/a in the fifth bullet point.

• Relative normalization between the transmit and the reflect component (T/R): We first discuss the results from the case: M_j → D_UXCL (⁠|$\rm j$| = RXT, MYT, CTOR, BOR). We find that T/R ≠ 1 (input). For the case M_BOR → D_UXCL, we get an exceptionally high value of T/R ≃ 6.38. Physically, the zeroth-order continuum across all models has the same shape for given a N_{H, los} (A_Fe = 1), whereas the scattered continuum is very different for different models. The tendency of T/R > 1 in CM-fits may be an attempt to get a better fit. Physically if the value of T/R is inconsistent with 1, it would imply variability of the coronal power law, however as demonstrated here this effect might be an artefact of model difference.

• Parameters of Torus Morphology:

  • Opening angle (⁠|$\rm cos \theta _{\rm o}$| or C_{frac, tor}): For M_BOR → D_CTOR we get cos θ_o ≃ 0.4 (θ_o ≃ 67^○) with corresponding θ_i ≃ 69.6^○. However in CTORUS the obscuring clumps are distributed in the region θ ≥ 60^○. This can be due to the difference in the total covering fraction of between the two torus. For the case M_BOR → D_UXCL in the coupled configuration (N_{H, los} = N_{H, eq} in this case, (Baloković et al. 2018)), we get a case where cos θ_o < cos θ_i which implies the line of sight does not intersect the main gas distribution of the torus, but returns a non-zero LOS absorption (⁠|$\rm \mathit{ N}_{H,los} \simeq 92$|⁠). This dichotomy can be misinterpreted as a different matter distribution (e.g. a stray clump, same N_{H, los} as the main gas distribution) lying along the LOS. In the uncoupled configuration of M_BOR → D_UXCL, the cos θ_o > cos θ_i, but the N_{H, los} ≃ 97 and N_{H, eq} ≃ 1300 heavily discrepant and becomes inconsistent with respect to the defined morphology. However, a physical interpretation for this case can be made, if the torus is clumpy, in which case |$\rm \cos \theta$| loses its meaning and hence can be interpreted as covering fraction, C_{frac, tor}.

  • c/a of RXTORUS: For M_RXT → D_MYT, c/a = 0.62 ± 0.03 (Fig. 11) which is discrepant with the fixed assumed value of 0.5 (effective θ_o = 60^○) in MYTORUS despite the fact that both are doughnut models. Given a set of parameters the scattered continuum of MYTORUS and RXTORUS have different shapes because of the radiative difference. The scattered continuum has more photon counts due to the presence of Rayleigh scattering in the E < 6 keV band and lesser photon counts at the CRH because of having E_cut fixed at 200 keV (Fig. 11). The parameter c/a thus adjusts itself to a higher value decreasing the effective θ_o to 51.7^○. Γ and θ_i were also reduced to 1.83 and 65^○, respectively (Fig. 11). For the M_RXT → D_UXCL case c/a ≃ 0.88, which can be considered consistent with the extended and random cloud distribution. This also led to θ_i for this case to take a comparatively low value, ∼37^○.

  • σ_o of UXCLUMPY: The value of σ_o follows different trends for data sets simulated with different models. For all the cases of M_UXCL → D_j, j = MYT, RXT, BOR, and CTOR, values of σ_o were systematically found to be ≤30^○ (Fig. 12) which means most of the clouds aggregate closer to the equator. For all the fits δ < 0.1, which is lower compared to the IM fits. In order to confirm that this trend is due to the effect of systematics (model difference) we performed an additional spectrum of MYTORUS spectrum for an exposure of t_XMM = 1.0 Ms and t_NuSTAR = 0.5 Ms and fit with UXCLUMPY (M_UXCL → D_MYT). We find that the median values and the errorbars (the diamond and circular markers in 12) are consistent with each other. This indicates that the trends for the high exposure simulations are similar to that of the typical exposure case of M_UXCL → D_MYT. This proves that the trends are systematic. We simulate data under UXCLUMPY using the best-fitting parameters obtained from M_UXCL → D_MYT and perform an IM fit. This is partially similar to the ‘surrogate’ method (Press et al. 2002). However, our analysis differs from Press et al. (2002), as we use a single ‘surrogate’ data set and our analysis of it is purely Bayesian. Consequently, the objective is comparison of the properties of the IM-posteriors obtained from the ‘surrogate’ data set with that of the CM-posteriors obtained from the real data set. We find that when the M_UXCL → D_MYT CM-fit returned δ ≃ 0.07 for σ_o, the IM fit for the ‘surrogate’ data set returned are wider with (Table 5) δ ≃ 0.46 for σ_o, which is ∼6 times higher.

  • Covering fraction (C_frac) of the inner ring of UXCLUMPY: The C_frac of the inner ring in the MCT regime returned posteriors with median values of C_frac < 0.1 (M_UXCL → D_j) or irregular posterior converging at the lower bound of the prior range implying C_frac ∼ 0 (Fig. 12) (all the other cases). We compare the posterior of C_frac from the CM-fit M_UXCL → D_MYT and the posterior obtained from IM-fit of the ‘surrogate’ data set. The posterior of C_frac in the CM-fit is irregular as it converges towards the lower edge (see Fig. 12). Similar to the σ_o the high exposure spectral fit (M_UXCL → D_MYT), returns trends consistent with the typical exposure spectrum, implying a systematic trend. The median is ∼0.01. However, the C_frac-posterior obtained from the same UXCLUMPY IM-fit of the ‘surrogate’ data set does not show irregularity. We find that C_{frac, median} ∼ 0.19 with δ ≃ 1 and R₉₀ ≃ 0.1 range is significantly large compared to that from M_UXCL → D_MYT.

  • N_cloud of CTORUS: The posteriors of the N_cloud for all the CM-fits, are irregular in the sense that the probability function has a minimum or is a monotonically increasing or decreasing function in the allowed prior range. The best values are either consistent with 2 or 10 or both (e.g. Fig. 10) all cases. In a real scenario, this can be an indication that CTORUS is a wrong model to apply in this case.

Table 5.

Parameter comparison for the cross-model fits (M_UXCL → D_j, j = MYT-MCT, CTOR-HCT, and BOR-HCT) and surrogate data sets simulated under the fitting model (UXCLUMPY) for the best-fitting parameters. θ_i and the parameters of morphology exhibit lower value R₉₀ compared to the intramodel-surrogate fit for UXCLUMPY.

Model parameters of MUXCL	MYT-MCT regime		CTOR-HCT regime		BOR-HCT regime
R90 values for CM-fit and surrogate data set
	DMYT	DUXCL (surrogate)	DCTOR	DUXCL(surrogate)	DBOR	DUXCL (surrogate)
σo(in °)	4.13	28.95	1.42	2.50	1.57	2.68
Cfrac of inner ring	0.04	0.26	0.004	0.027	0.008	0.034
θi	3.63	9.47	0.15	1.31	0.26	1.08

Table 6.

Table for input and the fit parameters for data simulated under torus and relativistic disc. Two cases are shown in this table. In the first case, only the torus (M_UXCL) is fit to the data where both the torus and disc components D_{UXCL + RELXILL} are present. In another case, the fit was performed with both the relativistic disc and the torus.

Parameters	AFe	log ξi	NH, los (cm−2)	Γ	σo (°)	Cfrac	θi (°)	χ2/dof
Input
MCT	1.0	1.0	100.0	1.9	45.0	0.4	60.0	–
HCT	1.0	1.0	500.0	1.9	45.0	0.4	60.0	–
MUXCL → DUXCL + RELXILL
MCT	–	–	102.0 ± 1.24	1.83 ± 0.03	|$32.0^{+8.1}_{-4.4}$|	|$0.57^{+0.02}_{-0.04}$|	|$46.6^{+5.6}_{-4.3}$|	0.978
HCT	–	–	|$547^{+58}_{-57}$|	1.72 ± 0.02	|$32.8^{8.9}_{-4.6}$|	|$0.47^{+0.01}_{-0.01}$|	|$24.0^{9.7}_{-14}$|	1.046
MUXCL + KD-XIL → DUXCL + RELXILL
MCT	|$0.61^{+0.14}_{-0.09}$|	|$0.62^{+0.64}_{-0.54}$|	97 ± 2	1.86 ± 0.03	|$48.6^{+11.5}_{-14}$|	|$0.34^{0.09}_{-0.1}$|	|$77^{+6.5}_{-14}$|	0.974
HCT	|$3.36^{+5.6}_{-2.41}$|	|$2.31^{+1.5}_{-2.0}$|	|$516^{+106}_{55}$|	|$1.9^{+0.04}_{-0.03}$|	|$51^{+12}_{-11}$|	|$0.44^{+0.06}_{-0.08}$|	|$50.6^{+11.4}_{-25}$|	1.035

4.2 Heavy Compton-thick regime

In the HCT regime we do not get good fits (compared to the MCT regime) for a lot of cases when models are applied in their consistent form. So in many cases we had to apply the fits with the zeroth-order continuum and the scattered continuum uncoupled and or ignoring the data where there is an emission line (one can alternatively add Gaussian lines), to get a good fit. We ignored the energy ranges connected with difference in the emission lines in the following specific cases: M_MYT → D_UXCL(1.16–1.31, 1.68–1.78, 3.5–3.61 keV), M_BOR → D_UXCL (2.15–3.26, 1.64–1.75 keV), M_RXT → D_UXCL (2.15–3.26, 1.64–1.75 keV), M_CTOR → D_UXCL (2.15–3.26, 1.64–1.75 keV), and M_MYT → D_RXT (1.6–1.9, 2.2–2.4 keV).

• Column density (N_{H, los}): The values of N_{H, los} were not consistently retrieved (e.g. Fig. 14b) for the CM-fits in the HCT regime. In the context of the input model, the strong extinction on the transmitted component worsens the recovery of N_{H, los} when wrong models are fit. For the HCT regime all models were simulated with N_{H, los} ≃ 500 and T/R = 1. The dissimilarity in the shapes of the CRH lead to bad fits in many cases when T/R = 1 (frozen) in the coupled configuration. So, T/R is kept free in most of the cases. We fit the doughnut (MYTORUS, RXTORUS) and biconical cut-out model (BORUS) to UXCLUMPY, in both coupled (N_{H, los} = f(N_{H, eq})) and the uncoupled configuration (where N_{H, los} ≠ f(N_{H, eq})). The doughnut fits were always significantly better in the uncoupled configuration but the N_{H, los} values are discrepant in all cases. The fits using BORUS seems practically uneffected in the uncoupled configuration. However, severe discrepancy and irregular posteriors in N_{H, los} are seen many cases (0.88 ≤ r_q ≤ 27.0, for the cases with regular N_{H, los} posteriors). There exist no specific trend methodology by which we can discern a correct value of N_{H, los}. The fitting is based on the attempt to adjust N_{H, los}, Γ and T/R such that it best mimics a CRH generated under a different morphology. This is unlike the MCT regime where the rollover and the FeK-edge from the zeroth order contunuum determines N_{H, los} and the effect of CRH on N_{H, los} relatively less.

• Photon-Index (Γ): Γ is discrepant (e.g. Fig. 14a) in most cases due to the varying shapes of the CRH. For all models the input value of Γ was 1.9. From the fits we find 0.67 < r_q < 21 holds for all tested model. The best value of r_q was obtained for r_q = 0.67 for M_RXT → D_UXCL (at Γ_out = 1.91) and the worse is obtained for M_CTOR → D_BOR (at Γ_out = 1.69). Thus in general it can be concluded that the recovered value of Γ is strongly model dependent.

• Relative normalization between the transmitted and the reflected component (T/R): When T/R was kept free, we find irregular posteriors and or severe discrepancy with the input value ([T/R]_input = 1 for all spectral simulations). One particular example case is M_RXT → D_UXCL where [T/R]_out = 4.92 in the uncoupled configuration. This results in stronger zeroth-order continuum in the fit (originally much weaker in the D_UXCL). The ∼ 5-fold increase in T/R effectively tries to fit the E > 10 keV region of the data. (Fig. 15). In other cases, T/R is reduced to very low values, thus arbitrarily decreasing the zeroth-order continuum (which is otherwise attenuated by N_{H, los}). This results relatively lower values of N_{H, los} (e.g. M_MYT → D_UXCL: [T/R]_out ∼ 0.5, N_{H, los} ≃ 200). Several other cases where inconsistency in T/R occur are: M_BOR → D_UXCL, M_CTOR → D_UXCL, M_CTOR → D_BOR etc. Physically, T/R ≠ 1 would imply variability of the coronal power-law; however as demonstrated here this can arise due to morphological difference resulting in varying shape of a CRH. The large discrepancy of T/R can result in heavy discrepancy in the intrinsic luminosity of the corona (L_corona).

• Parameters of morphology:

  • c/a of RXTORUS: For the case of M_RXT → D_UXCL(C_frac = 0) in the coupled configuration, c/a ≃ 0.11. The geometry was thus assumed to be more of an annulus than a doughnut but χ²/dof = 1.6. For the uncoupled configuration c/a returned values much higher, ∼0.71. However, the interpretation of c/a at face value is not possible. The geometrical interpretation of the results: the absorption happens through a Compton-thick absorber placed close to the axis of symmetry (θ_i ≃ 3^○) which is independent of a separate ∼ medium Compton-thick (N_{H, eq} ≃ 86, point 1) doughnut shaped gas distribution. For M_RXT → D_UXCL(C_frac = 0), c/a posteriors were bimodal, which resulted in bimodality in other posteriors. For both coupled and uncoupled configuration M_RXT → D_UXCL(C_frac = 0.4) returned c/a > 0.75 indicating an increased covering fraction of the torus. M_RXT → D_MYT also returned bimodal c/a-posterior.

  • cosθ_o or C_{Frac, tor}: In the case of M_BOR → D_CTOR we get cosθ_o ≃ 0.39, which is consistent with cosθ_i. The implications are similar to that of the fit obtained from the MCT regime. In the case of M_BOR → D_UXCL, the uncoupled configuration (N_{H, los} ≠ N_{H, eq}) returns a bimodal |$\rm cos \theta _{\rm o}$| or C_{Frac, tor} (one major and one minor peak in the posterior), which leads to bimodality in θ_i, A_Fe, and C_scpl. For this case we find cos θ_i < cos θ_o. Adding this up with the fact that N_{H, los} ≠ N_{H, eq} it is evident the otherwise ‘cos θ_o’-parameter makes more sense when interpreted as covering fraction C_{frac, tor} of a clumpy torus as described in the fifth point in Section 4.1.

  • σ_o of UXCLUMPY: M_UXCL → D_CTOR returns σ_o ≃ 14^○ with δ ≃ 0.1 and M_UXCL → D_BOR returns σ_o ≃ 12^○ and δ ≃ 0.1. In both the cases the joint analysis of XMM–Newton and NuSTAR returned fits for which |$1.2 \le \, \chi ^2$|/dof ≤1.3 (mainly due to discrepancy in the emission lines). We simulate ‘surrogate’ data sets under the best-fitting parameters of UXCLUMPY under both M_UXCL → D_CTOR and M_UXCL → D_BOR best fits parameters and carry out IM-fits. We find that the IM-fits in the surrogate data sets return wider posteriors (Table 5). For σ_o, δ_IM/δ_CM is 1.75 for BORUS and 1.77 for CTORUS. Both the fits M_UXCL → D_RXT and M_UXCL → D_RXT returned very low values of σ = 7^○.

  • C_frac of the inner ring of UXCLUMPY: Unlike the MCT regime, values of C_frac returned from the M_UXCL → D_CTOR and M_UXCL → D_BOR were clearly non-zero and assumed significantly higher values (C_frac ≃ 0.3 for both) with tighter constraints δ ∼ 0.01 for M_UXCL → D_BOR and δ ∼ 0.001 for M_UXCL → D_CTOR. For the NuSTAR only fits as described in the previous point, the values of C_frac and δ are similar to those of the joint fits for M_UXCL → D_BOR. For the M_UXCL → D_CTOR fits the C_frac value decreased but the δ increased to 0.5. The ‘surrogate’ data sets under IM fits returned broader posteriors (Table 5). For C_frac, δ_IM/δ_CM is 4.04 for BORUS and 6.67 for CTORUS. For both the RXT and MYT fits the posteriors are similar in shape and the median values are similar.

  • N_cloud of CTORUS: We carried out the CM fits M_CTOR → D_j where j = MYT, UXCL, and BOR. For M_CTOR → D_UXCL and D_MYT we do not get a good fit (⁠|$\rm 1.2 \le \chi ^2/dof \le 2.0$|⁠) mostly because of the difference in the emission line profiles between UXCLUMPY and CTORUS. We get an irregular posterior for N_cloud converging towards the higher extreme edge of the prior range for both cases. For M_CTOR → D_BOR, the N_cloud posteriors are extremely narrow.

Figure 14.

Cross- and intramodel fits for the case M_j → D_UXCL, j = MYT, RXT, CTOR, and BOR models. The labels: -c and -uc refers to coupled and uncoupled configuration, respectively. The discrepancies in Γ are due to both flattening and steepening of Γ. The discrepancies in N_{H, los} are also large amounting to a maximum of |$\sim 100 {{\ \rm per\ cent}}$|⁠.

Figure 15.

M_RXT → D_UXCL fit for two cases: Z–S coupled with χ²/dof ≃ 1.4 and log BF ≃ −364.99; Z–S un-coupled: with χ²/dof ≃ 1.16 and log BF ≃ −148.56. The Bayes Factors are calculated using the log |$\rm \mathit{ Z}_{\rm UXCL}$| value obtained from the UXCLUMPY intramodel fit. Most of the discrepancy in the spectra happens in the continuum.

5 RESULTS: DETECTABILITY OF AN ADDITIONAL RELATIVISTIC DISC-REFLECTION COMPONENT

In a Compton-thick AGN, the transmitted and reprocessed emission of the torus is the dominant part of the spectrum. However, there may be an additional component present in the spectrum originating due to relativistic reflection from the inner accretion disc. The question thus is, would it be possible to detect such a component from spectral model fits? To answer the question, we simulate a data set that has both the CRH-dominated torus component and a relativistic disc-reflection component. We use UXCLUMPY to simulate the torus and RELXILL (García et al. 2014) to simulate the disc-reflection component (D_{UXCL + RELXILL}). Since the disc reflection component originates close to the ISCO, it should be heavily absorbed by the torus column. Thus the model expression in xspec terminology is:

MODEL = TBABS×(APEC(1)+APEC(2)+ UXCLUMPY-SCPL + ZTBABS×CABS×RELXILL + UXCLUMPY)

The ZTBABS and CABS N_{H, los} was tied to UXCLUMPY as the obscurers are one or multiple torus clouds. The parameters of RELXILL like Γ and θ_i were tied to the UXCLUMPY as the torus and the disc are illuminated by the same coronal power law. In the flux range of 4–50 keV, the ratio of the flux of the disc absorbed through the torus column with respect to the total torus emission F_{4–50, Disc}/(F_{4–50, trans} + F_{4–50, refl})_Torus is 0.36–0.38, which makes the total torus emission (comprising of the torus CRH and the transmitted component) the dominant feature of the spectra. In the published UXCLUMPY model, the CRH from the torus is calculated assuming a power law. But when a disc reflection is present, the torus should reprocess both the power-law and disc-reflection component. It was not possible for us to implement this in the current UXCLUMPY setup. This will put some limitations on the accuracy of our test but we do not expect it to significantly change the results. When our fitting model has only torus component (UXCLUMPY) and the data has both the disc RELXILL and torus component the notation would be M_UXCL → D_{UXCL + RELXILL} and in the converse case, the disc is modelled with KDBLUR(XILLVER), the notation is M_{UXCL + KD-XIL} → D_UXCL. During analysis, the prior of the relative normalization of the disc-component C_KD-XILL was set to be a log-uniform prior extending from 10⁻¹⁰ to 1. All other parameters of the disc that were kept free (ionization parameter logξ_i and abundance A_Fe) were given uniform priors extending accross the entire prior range. The methodology of analysis of the data is described in the following bullet points:

• Torus component only: We fit the data with torus component only i.e. M_UXCL → D_{UXCL + RELXILL} to see if it is possible for only a torus modeled by UXCLUMPY to replicate the spectra. We keep the APEC components and E_cut frozen to their input values. We find that the spectral shape was replicated and there no significant model or data excess in the residuals for both MCT and HCT regime. The overplot of the input-model and fit model for the HCT regime is shown in Fig. 16. At the median values of the parameter posteriors, we obtain χ²/dof = 1.046 for HCT and 0.978 for MCT. The parameters like Γ, C_frac of UXCLUMPY and θ_i showed large discrepancies with the input. σ_o was consistent with the input. The fit returned a lower value of Γ and a higher value of C_frac (see Table 6).

• Torus and disc component: We fit the simulated data with UXCLUMPY and an attenuated disc component modelled by ZTBABS × CABS × kdblur⁶ (XILLVER⁷). This produced a better fit compared to the previous case, with χ²/dof ≃ 1.035 for HCT and 0.974 for MCT, with no significant residuals. The XILLVER parameters like A_Fe and χ_i show very broad posteriors, which is attributable to the heavy extinction through the torus column.

Detecting an additional relativistic disc reflection component which is heavily attenuated in the torus absorbing material with certainty is impossible with simple χ²/dof values. In all the fits with only the torus model (UXCLUMPY), the model replicates the data with almost indistinguishable values of χ²/dof. But the low Bayes factor (see Section 6 for definition), BF = |$\rm \mathit{ Z}_{torus}/\mathit{ Z}_{torus+disc} \simeq 1.21 \times 10^{-19}$| for HCT and 9.25 × 10⁻⁸ for MCT clearly hints at the existence of a relativistic disc reflection component. We also test the possibility for false detection of a disc component in the HCT regime and carry out the M_{UXCL + KD-XIL} → D_UXCL fit. We found that the parameters of UXCLUMPY got recovered, but the parameters of XILLVER showed irregular posteriors. Specifically, the tendency of the C_KD-XILL posterior-distribution to converge towards the lower limit (∼10⁻⁸) indicated the absence of the disc reflection component. The log-evidence values for the M_{UXCL + KD-XIL} → D_UXCL and M_UXCL → D_UXCL are log|$\rm \mathit{ Z}_{torus+disc} \simeq -652.19$| and log|$\rm \mathit{ Z}_{torus} \simeq -653.67$|⁠, respectively. Thus, the bayes factor is BF = |$\rm \mathit{ Z}_{torus}/\mathit{ Z}_{torus+disc} \simeq 0.033$|⁠, the value being greater than our threshold of 0.01 (Section 6) for model distinction. Nevertheless, the results of M_{UXCL + KD-XIL} → D_UXCL is consistent with the absence of a relativistic disc refelction.

Figure 16.

Cross model fits: overplot of UXCLUMPY input model and fitting models (a) in the MCT regime (b) for the HCT regime with the Compton-thick ring absent (C_frac = 0) (c) Overplot of CTORUS input and M_UXCL → D_CTOR best-fitting model. The zeroth-order continuum adjusts itself to replicate the CRH of CTORUS (d) Overplot of input model comprised on a weak disc component(RELXILL) and torus component(UXCLUMPY) and fitting model with only a torus model(UXCLUMPY).

6 RESULTS: DEPENDENCE OF EVIDENCE ON FLUX LEVELS

To study the effect of flux levels on model evidence, we simulate 10 spectra (S_n, n runs from 1 to 10) for 10 different flux levels under UXCLUMPY and CTORUS for an intrinsic 2–10 keV flux of 0.83/n mCrb in the HCT regime. The exposures on all the simulated spectra are 100 ks on XMM–Newton and 50 ks on NuSTAR. All spectra are grouped into 30 counts bin⁻¹. We run BXA sims on the simulated spectra for both: (A) joint XMM–Newton and NuSTAR (B)NuSTAR only cases. Our Bayesian fitting strategy can be summarized in the following representation: M_j → D_UXCL and M_j → D_CTOR, where j = UXCL, CTOR, and BOR. In our work, Bayes factor is defined as the ratio of the Bayesian evidence of a fitting model (Z_{fitting model}) to that obtained from an intramodel fit (Z_input): |$\rm BF = Z_{\rm fitting-model}/Z_{input}$|⁠. We plot the Bayes factor values to the corresponding flux levels (Fig. 17). For effective distinguishing of models, we set the following condition on Bayes factor (BF): |$\rm BF \lt 10^{-2}$|⁠, which practically means that the flux value below which the |$\rm log BF$| value rises above −2 for the first time is considered as the lower limit on flux, below which effective distinguishing of correct from the wrong model is not possible. We expect that joint fits are more effective in distinguishing models. Our results confirm this from the fact that we can distinguish the models more effectively down to an intrinsic flux value of ∼0.1 mCrb for joint XMM–Newton plus NuSTAR fits. For the NuSTAR only fits, for which this lower threshold value was slightly higher than ∼0.1 mCrb (Fig. 17b and d). It is also worth noticing that systematically lower values of logBF and local fluctuations in logBF (e.g. as seen in the NuSTAR-only fits; Fig. 17) can push the threshold to higher limits (grey vertical line at flux level f/3 in Fig. 17d). These fluctuations which might have occurred during the data simulation or detection process are more likely to influence the evidence values at the lower flux regimes or lower bandwidth data, thus reducing the precision in evidence-based model distinction. Thus, it is safer to assume more conservative thresholds determined from synthetic data simulations, when applying Bayes factor for model distinction in case of real data. We discuss these aspects in detail in Section 7.

Figure 17.

log Bayes Factor versus flux level plot: Spectra simulated for XMM–Newton and NuSTAR exposures of 100 and 50 ks, respectively, with an intrinsic flux level of f_n = 0.83/n mCrb. (a) and (c) demonstrate the variation of evidence values with flux when both XMM–Newton and NuSTAR data are used. (b) and (d) demonstrate the same, but considering NuSTAR only data. The blue and magenta dashed lines denote the flux level at which the relation |$\rm log \mathit{ Z}_{\rm input} - log \mathit{ Z}_{\rm fit} \rm \ge 2$| holds for the wrong model. The black and green dashed lines marks the logBF = 2 and logBF = 0 level, respectively.

7 DISCUSSIONS

7.1 Guidelines to the X-ray community

7.1.1 Parameter estimates

We perform joint analysis of synthetic data simulated under XMM–Newton and NuSTAR. The baseline flux level for all the input models was ∼0.5 mCrb in the 2–10 keV band. The exposure times on XMM–Newton and NuSTAR were 100 and 50 ks, respectively. For XMM–Newton EPIC-pn, given the flux and exposure time in the 2–10 keV band, the average photon-count is ∼8.2 × 10⁴. From the results of our analysis of the synthetic data, we assert the following prescriptions for the parameter posteriors:

• From the cross-model fits, it is quite clear that the photon-index (Γ) cannot be determined with certainty when a wrong model is applied. Almost all the posteriors are regular-monomodal distributions with very small statistical errors, but the discrepancy (represented by the r_q values) between the input and fit values are extraordinarily high (ΔΓ ≃ 0.3 and r_q as high as 16). This discrepancy in Γ for different models is also reported by Buchner et al. (2021). Even multiple models returning similar values of Γ should not be taken at face value (as discussed in Section 4). Our results indicate that precision in the value Γ is no guarantee for its correctness. In such a scenario, where the systematic discrepancies originating due to model-difference cannot be corrected, we recommend that users apply caution in the interpretation of Γ (e.g. when used to estimate intrinsic luminosity or accretion rate relative to Eddington), particularly if very low values of Γ are obtained. Alternate, indirect methods to determine the value of Γ may be preferred. For example, Brightman et al. (2013) noted an empirical relation between Γ and FWHM of optical broad emission lines in X-ray selected, broad-line radio-quiet AGN. However, such a sample is biased to be X-ray-unobscured, and this relationship might not hold in X-ray obscured AGNs. The community would thus benefit from exploration into whether any empirical relation between Γ and other measurable properties of Compton-thick AGN exists.

• Consistent values of LOS column density (N_{H, los}) were returned in the medium Compton thick (MCT) regime, because of the strong zeroth-order continuum. In the heavy Compton thick (HCT) regime, the values of N_{H, los} were recovered but were less constrained as seen from most of the intramodel fits (IM-fits). However, when wrong models were applied (CM-fits) many cases returned discrepant values and irregular posteriors. The determination of the precise value of N_{H, los} starts to become difficult in the HCT regime, because the suppression of a zeroth-order continuum by the scattered continuum. However also in the HCT regime, in the CM-fits there are cases when the intensity of the zeroth-order continuum is modulated by the T/R-parameter (arbitrarily taking values > or < 1) to adjust the shape of the CRH (Fig. 15). The large discrepancy of T/R can result in an equally sized mis-estimation of the intrinsic accretion luminosity (L_corona). (L_{X, corona}), which can be severely discrepant (L_{corona, fit} upto ∼ 10 times L_{corona, input} in our simulations).

• The IM-fits returned comparatively narrower posteriors (Fig. 8) in the HCT regime for most parameters of torus morphology. The constraints on the parameters of morphology are drawn from the E < 6 keV tail, Fe K edge, and CRH of the scattered continuum which are dominant in the HCT regime. However, these comparative levels of improvements in the constraints in the HCT regime w.r.t. to MCT can be dependent on the exact spectral shape. However, it can be inferred that under the correct model assumption parameters of morphology on average can be better constrained in the HCT regime, as the HCT regime has lower ‘contamination’ from the zeroth-order continuum. However, the behaviour of the posteriors of the morphological parameters shows different trends for different models for intramodel fits and cross-model and are open to interpretation:

  • Continuous torus: In the continuous torus model there are only two morphological parameters: c/a for doughnuts and θ_o or cos θ_o or C_{frac, tor} in the biconical cutouts. The doughnut thickness parameter c/a (of RXTORUS) shows regular posteriors (Section 2.4 for definition) in all cases. However, low values of c/a (e.g. c/a < 0.2) were observed for the baseline (coupled) configuration, which indicates an annular geometry rather than a doughnut. Meanwhile for the uncoupled configuration one interpretation of the value of c/a can indicate the reflector geometry generating the scattered component (see Section 4.2). Several other interpretations with different setups for an otherwise doughnut model are explained in Yaqoob (2012). cos θ_o or C_{frac, tor} of BORUS returns posteriors which are regular for most cases. However for some fits to wrong models in the coupled configuration (N_{H, los} = N_{H, eq}) we get θ_o > θ_i. This suggests an unobscured configuration. However these cases returned a N_{H, los} consistent with a Compton thick obscured case. For our case, we know that the applied model is incorrect. However in real situation where the input spectrum is unknown and these class of solutions are encountered there may a create potential chance of misinterpretation of an otherwise wrong solution (see the point on cos θ or C_{frac, tor} in the Section 4.1).

  • Clumpy torus: In the clumpy torus models there are three morphological parameters: σ_o and C_frac of inner ring of UXCLUMPY and N_cloud of CTORUS. When UXCLUMPY is the wrong model (CM-fits) σ_o systematically showed narrower posteriors (lower values of R₉₀ and δ) compared to the case when it is the correct model (IM-fits) at the same flux levels. The covering fraction (C_frac) of the inner C-thick ring of UXCLUMPY returned irregular posteriors (Section 2.4 for definition) with C_frac ∼ 0 for majority of the tested cases in the MCT regime. However in the HCT regime the posteriors of C_frac were monomodal with comparatively higher values (C_frac > 0.1). This could be interpreted as the tendency of the UXCLUMPY geometry to adjust the overall covering fraction of the clump distribution to match the shape of the data generated under another morphology. N_cloud of CTORUS returned mostly irregular posteriors when CTORUS is the wrong fitting model, which had either a minimum or converged at one of the extreme edges of the priors. For the intramodel fits, however, the N_cloud posteriors are wide but mostly regular.

  Although it might not be a generic trend there might exist differences in the nature of the posteriors for a wrong and correct model. Thus in real data analysis a rigorous test of the behaviour of the parameters of the fitting model becomes a necessity. This can be done by simulating synthetic data under the fitting model for the values of the best-fitting parameters and then performing an intramodel fit. The posteriors properties (e.g. δ, R₉₀) for the real data and the intramodel fit can then be compared.

• Posteriors of θ_i are generally regular monomodal in case of application of the correct model (IM fits). However, in several cases the proportional errors are large. In general, the constraints on θ_i are dependent on the constraints on other morphological parameters, which are better determined in the HCT regime with a dominant scattered continuum allowing better constraints. However, when the applied model is wrong (CM-fits), in several cases we got irregular posteriors (converging at the edge) especially those involving the fitting of UXCLUMPY. In general, θ_i remains a parameter dependent on morphological assumptions and hence should not always be accepted at face value.

• As described before a significant fraction of cases where a wrong model is applied (CM fits), returned irregular posteriors for some of the parameters [e.g. θ_i, N_cloud, C_frac etc. (see Table 1) for a summary of the parameters]. However, for a significant number of fitting models, the best value (typically the median) derived from the posteriors of the cross-model fits returned χ²/dof < 1.2, which is otherwise acceptable. Also in most of these cases, it is possible to improve the value of χ²/dof significantly by adding emission lines or by decoupling the zeroth-order and the scattered continuum (see Section 4.2), etc. This might return meaningful values of certain parameters (e.g. N_{H, los}), but the morphological and radiative aspects remain open to (mis-)interpretation (e.g. point-3 of Sections 4.1 and 4.2). The danger of misinterpretation is less in the case for multinest as in addition to χ²/dof estimates we can apply Bayes-Factor to weigh up the model preferences and explore the nature of intramodel posteriors.

We thus strongly recommend that the users perform their data simulations and use global Bayesian algorithms for data analysis. The methodology used here and the exercise undertaken in this work can be treated as an example. We advocate simulating data under selected models of choice to understand potential limitations both in cases of preparing observing proposals and when modeling already obtained data. We leave it to the readers to use the responses/ARFs of the selected instrument(s) as appropriate, and perform intramodel and/or cross-model fits as needed to understand the behaviour of the posteriors for the given models, instrument, exposure time, etc.

7.1.2 Limitations of using Bayes Factor values

From the results in Section 6, it is clear that values of log-evidence (logZ) and log-Bayes factor (log BF) are strongly affected by the number of degrees of freedom, which is a strong function of the number of data bins, and availability (or lack thereof) of certain instruments and bandpasses. For a given flux/exposure time, maximizing the energy bandpass used by, for example, performing joint NuSTAR and XMM–Newton fits instead of fitting just one instrument maximizes the ability to distinguish between models based on Bayes Factor. In the case of joint NuSTAR and XMM–Newton fits to the synthetic Compton-thick AGN spectra (simulated under the selected models; Section 6), adopting a threshold log-Bayes factor of ≳2 was effective in model distinction for cases with 2–10 keV fluxes brighter than ∼0.08 mCrb (Fig. 17), for 50/100 ks exposures with NuSTAR/XMM–Newton, respectively. For the NuSTAR-only fits the availability of fewer channels (data points) brings the logBF (Bayes factor) value closer to 0 (1), for a given intrinsic flux value. Additionally, as demonstrated in Section 6, ‘local’ random fluctuations in values of logBF occur (Fig. 17), but whether these fluctuations will impact distinction of a wrong model from a correct model depends flux or availability of data in a certain band. For our examples (Section 6 and Fig. 17), for the joint fits, the values of logBF are more negative and lie comfortably below the selected threshold (up to a flux range of f/9 or 0.09 mCrb) despite random fluctuations. Meanwhile, for the NuSTAR-only fits the values of logBF are closer to the threshold and thus run more into the risk of crossing into the logBF > 0 zone under the influence of random fluctuations leading to a false inference on the morphology. The grey and vertical dashed line in Fig. 17(d) marks the flux value where logBF for M_UXCL → D_CTOR comes dangerously close to the logBF-threshold (but does not cross it) for a NuSTAR-only fit. Thus, there exists a possibility that for a given threshold on logBF, the threshold on flux required for an effective model distinction can get pushed to a higher value in the absence of a certain band/instrument and or presence of random fluctuations in values of logBF. Hence, these threshold values are approximate estimates and should not be considered as absolute. It should also be noted that these results from our examples are based on the assumption that the soft band emission (e.g. emission modelled with APEC) is modelled perfectly, which might not be the case in real situations(see Section 7.2). This is because of the involvement of several complicated factors related to measurement, contamination from other sources in the field of view and physical modelling of the system, we advocate that the readers perform their simulations setting a reasonable limiting value of logBF, to estimate the appropriate threshold on intrinsic flux necessary for model distinction.

7.2 Limitations of our analysis

• The spectral shape can determine counts in a band and various input parameters used for data simulation. For example, different values of N_{H, los}, different levels of contribution from the scattered power law, etc. can contribute to enhancement or suppression of torus features in the 2–10 keV band. This might also affect the relative spectral counts collected by XMM–Newton and NuSTAR given an exposure time, thus resulting in variation in the constraints. Effects of variation in spectral shape on parameter constraints were shown for a few parameters like N_cloud (Fig. 6). Comprehensive testing of all possible cases of spectral shapes for all models is out of the scope of this work. However, in the case of real data analysis, these spectral shapes should be taken into consideration.

• In real data reduction, the given spectrum has unknown parameters and we model it as best we can. But there will always be some discrepancy between the modelled parameters and the real spectral parameters. This consequently results in statistical and systematic variations in the location estimates (e.g. median) of the parameter posteriors during the data generation process. A similar situation is encountered for simulated data. To undertake such a study one can simulate a very large number of spectra (of order 10⁴) (e.g. González-Martín et al. 2019) using fakeit in xspec (e.g. Markowitz, Reeves & Braito 2006) and then carry out Bayesian analysis in bulk. We performed such a study on a smaller scale in Section 3.3, wherein we simulated and fit 100 spectra. Given the heavy computational requirements for the study, we limit ourselves to only 100 spectral fits. The distributions of the posterior properties can be studied better with a much larger number of simulated data sets.

• The models used here differ in two aspects namely, the assumed morphology and the assumptions on radiative physics. The results of our simulations have been limited by the simultaneous occurrence of both radiative and morphological differences between the models. Thus, we could not effectively differentiate between the effects of the radiative and the morphological differences in our current study. To test solely the effect of the morphological difference of the torus, it is necessary to set the radiative properties across different morphologies to be constant, so that the various properties of the CRH are just the effect of torus morphology only. Therefore, only one radiative transfer code e.g. XARS should be used to simulate different spectra for different morphology which in turn should be used for model testing. Conversely, to test the variation that radiative physics, different radiative transfer codes like XARS, REFLEX, GEANT, etc. could be used to simulate the spectra for a fixed morphology.

• In our analysis, we assumed only simplistic situations while modelling the emitted radiation from a Compton-thick obscured AGNs. More specifically we did not take into account any additional emission from other X-ray bright sources which might be present in the host galaxy. In a real scenario, there might exist significant contamination from other X-ray bright sources like ULXs, X-ray binaries, etc. This has been demonstrated in Arévalo et al. (2014) for the case of the Circinus galaxy, where several X-ray bright sources contaminate different energy bands at different levels. In our work, we have demonstrated the energy/band dependence of certain parameters. Thus, contamination can lead to wrong estimates of the parameter. It is thus required that the contaminants are accounted for and are modelled correctly.

• In the future, it is expected that more physically self-consistent torus models and complex-torus models (shown schematically in e.g. Esparza Arredondo et al. 2021) will be developed, incorporating additional morphologies, cloud distributions, and physical processes. In addition, the community can look forward to data from new X-ray missions, including XRISM and ATHENA: each of these missions will feature a calorimeter that can resolve subtle features in fluorescent emission line complexes such as the Compton shoulder (e.g. Yaqoob & Murphy 2011; Tanimoto et al. 2019). Even presently, data from currently active missions other than XMM–Newton and NuSTAR can be used for analysis. Thus, performing additional tests similar to that performed in this work is a necessity, when one uses a model and data from an instrument not been tested here.

8 SUMMARY

In this project we test the reliability of the estimates of the model parameters for an X-ray obscured Compton-thick AGNs. We simulate synthetic data under six models with the instrument functions of XMM–Newton and NuSTAR at a flux level of ∼0.5 mCrb in the 2–10 keV band. The exposure times were 100 ks on XMM–Newton and 50 ks on NuSTAR. Using Bayesian analysis techniques we perform intramodel (IM) and cross-model fits on the synthetic data using these models. We investigate the level of degeneracy between the parameters, the ability to distinguish models and the different ways a wrong model fitting results in incorrect interpretation of the morphological and central engine parameters (e.g. photon index). We point out the challenges and limitations concerning the analysis of single epoch spectra.

In the intramodel fits, we find that the parameters are mostly recovered. However many morphological parameters showed very wide posteriors covering more than |$30{{\ \rm per\ cent}}$| of the prior range in the medium Compton-thick regime (N_{H, los} < 2 × 10²⁴ cm⁻²). Constraints on N_{H, los} worsen and the constraints on the morphological parameters improve in the heavy Compton-thick regime (N_{H, los} > 2 × 10²⁴ cm⁻²). We also demonstrated the effect of a particular energy band (or instrument) which is necessary for determining reasonably good constraints on the parameters.

We also discuss a few ways in which a wrong model can fit a given data set and lead to false conclusions on the morphology along with other properties (e.g. N_{H, los}/N_{H, eq}) of the obscurer and or reflector. When a wrong model is fit to a particular data set, the photon index (Γ) of the coronal power law (a common parameter across all models) in most of the cases returned posteriors showing that Γ is constrained with a proportional error of <3 per cent but with significant discrepancy with the input. However for wrong model fits in case of the parameters which are present in the fitting model but not in the model used for data simulation, a different method can be adopted. For a few cases involving fitting with UXCLUMPY viz. M_UXCL → D_j where j = MYTORUS, CTORUS, and BORUS, we compare the posteriors obtained from the CM-fit to an that obtained from and IM-fit to a ‘surrogate’ data set (simulated under fitting model with the median values obtained from the CM-fits). The analysis of the ‘surrogate’ data set returned broader posteriors compared to the CM-fits for the morphological parameters σ_o and C_frac of the inner ring (see model description). It is thus established that the precision in the value obtained from a fit is no guarantee for its correctness and does not necessarily indicate the application of a model consistent with the real morphology. Our results from the simulations thus indicate the dangers of accepting fit results at face value.

The parameter spaces of Compton-thick AGN with limited data are complicated and sometimes multimodal, and thus demand global parameter exploration algorithms. With the use of Bayes-factor from bxa/multinest runs it is possible to rule out the least preferred models for a given data set, at least at the higher end the available data quality. However, the use of synthetic data (simulated under the models in use) to assess the possibility of random fluctuations in the values of Bayes factor, which can potentially lead to false model conclusions (Section 6, 7.1.2, and Fig. 17), is recommended given an instrument combination and or flux regime. We thus conclude that measuring the properties of Compton-thick obscurers or reflectors with single epoch data obtained jointly from XMM–Newton and NuSTAR bands is not straightforward, even without the complicated effects of contamination. However, a study with multiepoch data can help us study variability (to detect variations in N_{H, los}) and thus partial conclusion on the morphology can be drawn (e.g. Markowitz et al. 2014; Baloković et al. 2018; Laha et al. 2020; Esparza Arredondo et al. 2021). Thus time-variability analysis coupled with spectral analysis can yield better constraints on morphology.

The results of all the analysis carried out in this work thus suggest the importance of performing synthetic data simulation and analysis alongside real data using the instruments and the models of interest. This will help understand the potential limitations of the analysis, the models and the instruments. Hence, it is recommended that the users of these models perform similar simulations considering this work as example and noting corresponding potential limitations.

SUPPORTING INFORMATION

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

The authors thank Prof. Dr. hab. Agata Różańska for her valuable feedback on the scientific results and their broader significance in the X-ray astronomy. The authors also thank the anomymous referee for the insightful comments. This research has made use of data and software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC. This work was financially supported by Polish National Science Center (NCN) grant number: 2018/31/G/ST9/03224. AM also acknowledges partial support from NCN grant number: 2016/23/B/ST9/03123.

DATA AVAILABILITY

The corner-plots of all the analysed cases can be found in https://users.camk.edu.pl/tathagata/tsaha1996.html. Any other data products will be made available on reasonable request.

Footnotes

1

https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/node234.html

2

http://mytorus.com/mytorus-instructions.html

3

https://heasarc.gsfc.nasa.gov/xstar/docs/html/node106.html

4

http://atomdb.org/

5

The 2D posterior figures are produced using the package corner (Foreman-Mackey 2016). All other plots are made using the package matplotlib (Hunter 2007).

6

Laor (1991).

7

García et al. (2013).