© The Authors 2023

1 Introduction

Large amounts of gas and dust can be found in the central regions of most active galaxies. These ambient media play a crucial role as they provide the accretion reservoir powering active galactic nuclei (AGN), and they reprocess the optical, ultraviolet (UV), and X-ray emission of the central engine. Through these cir-cumnuclear media, the large-scale AGN feedback on the host galaxy can be studied, and the conditions close to the super-massive black hole can be probed (Ramos Almeida & Ricci 2017). Yet, the detailed characteristics of these regions and the physics behind them remain poorly understood. According to the AGN unification scheme, dusty gas exists as a toroidal structure in the equatorial plane, which explains the observed variety of AGN types as a distribution of system orientations relative to the observer (Antonucci 1993; Urry & Padovani 1995; Netzer 2015). This putative dusty torus is assumed to absorb 40–70% of the accretion disk emission, which is re-emitted at thermal infrared (IR) wavelengths peaking at about 10 to 20 µm (see e.g. Netzer 2015; Stalevski et al. 2016, and references therein). The successful modelling of observational IR spectra using dusty torus models supports this unified picture of obscuration in local Seyfert-type AGN (Siebenmorgen et al. 2005; Hao et al. 2005; Spoon et al. 2007; Schweitzer et al. 2008; Mendoza-Castrejón et al. 2015; Hatziminaoglou et al. 2015).

However, recent high angular resolution mid-IR observations of AGN suggest that the distribution of obscuring dust is also extended in the polar direction, as opposed to the classical dusty torus paradigm (Hönig et al. 2012; López-Gonzaga et al. 2016; Asmus et al. 2016). These polar dust features are now firmly established in two nearby active galaxies, NGC 1068 (Jaffe et al. 2004; Wittkowski et al. 2004; Raban et al. 2009; López-Gonzaga et al. 2014) and the Circinus Galaxy (Tristram et al. 2007, 2014; Stalevski et al. 2019; Isbell et al. 2022), while similar mid-iR extensions have been confirmed for other AGN (Hönig et al. 2013; Asmus 2019; Leftley et al. 2019). The physical origin of these polar structures could be radiation-driven winds that are launched close to the dust sublimation radius (Hönig et al. 2012), as demonstrated in theoretical calculations by Gallagher et al. (2015), Chan & Krolik (2016), Wada et al. (2016), and Vollmer et al. (2018). Furthermore, detailed radiative transfer modelling suggests that the AGN-obscuring medium has a complex structure with clumps and filaments, possibly forming an intricate two-phase medium of clump and inter-clump regions (Stalevski et al. 2012, 2016).

A promising way to study gas and dust in active galaxies is to look at their X-ray emission, as the interactions between X-rays and the circumnuclear medium produce characteristic spectral features that contain information on this obscuring body of gas and dust (Ricci et al. 2014; Ichikawa et al. 2019). Observations reveal strong X-ray emission emerging from most AGN, which is likely produced in a corona of hot electrons close to the central black hole, where optical and UV disk photons are Compton up-scattered to X-ray energies (Haardt & Maraschi 1991). This X-ray component is detected even for heavily obscured sources, which opens an interesting observational window on AGN obscuration (Ricci et al. 2017). While the transmitted X-ray continuum displays absorption features linked to the intervening gas density and abundances, most information is encoded in the reflected X-ray emission, which is now used to constrain the putative torus geometry and its physical properties (Gupta et al. 2021; Yamada et al. 2021; Osorio-Clavijo et al. 2022). The main signature of X-ray reflection in obscured AGN is the appearance of various fluorescent emission lines, the most prominent one being the narrow Fe K_a line at 6.4 keV (Mushotzky et al. 1978; Nandra et al. 1989; Pounds et al. 1989; Fukazawa et al. 2011). In addition, X-ray scattering produces characteristic spectral features such as a Compton reflection hump peaking at about 30 keV, and a Compton shoulder at the low energy side of bright fluorescent lines, which are believed to form powerful probes on the circumnu-clear material in AGN (Matt 2002; Matt et al. 2003; Yaqoob & Murphy 2011; Odaka et al. 2016; Hikitani et al. 2018).

In order to fully exploit the diagnostic power of these reprocessed X-ray features, one has to model the relevant X-ray absorption, scattering, and re-emission processes in great detail. Furthermore, to make realistic model predictions for AGN circumnuclear media, one needs to perform radiative transfer simulations in three dimensions. The detailed 3D substructure of the obscuring medium has been shown to have a strong effect on the dust modelling results in the IR-to-UV wavelength range, which was investigated with radiative transfer simulations based on smooth torus models (Fritz et al. 2006), clumpy torus models (Nenkova et al. 2008a,b; Hönig et al. 2006; Schartmann et al. 2008), two-phase filamentary torus models (Stalevski et al. 2012), and most recently two-component disk-plus-wind models (Hönig & Kishimoto 2017; Stalevski et al. 2017, 2019; Hönig 2019). Yet, the X-ray signatures of these 3D torus models have not been studied to the same level of sophistication as it has been done in the IR. The introduction of advanced X-ray torus models consistent with the 3D IR dust models would allow for a joint analysis of the X-ray and IR bands, which could break some degeneracies that are encountered when dissecting dusty structures using IR observations only (Farrah et al. 2016; Ogawa et al. 2019; Tanimoto et al. 2019; Esparza-Arredondo et al. 2021). Finally, the advent of X-ray microcalorimeter observations with XRISM/Resolve (Tashiro et al. 2020) and Athena/X-IFU (Barret et al. 2018) could revolutionise our understanding of AGN obscuration, given that detailed X-ray models are available that can be used to leverage the high energy resolution to constrain their precise distribution of gas and dust.

However, reprocessed X-ray spectra are often fitted with 1D XSPEC models (Arnaud 1996) such as PEXRAV (Magdziarz & Zdziarski 1995) and PEXMON (Nandra et al. 2007), which cannot capture the 3D complexity of the circumnuclear medium. More advanced X-ray spectral models are available (Murphy & Yaqoob 2009; Ikeda et al. 2009; Brightman & Nandra 2011; Odaka et al. 2011, 2016; Liu & Li 2014; Furui et al. 2016; Paltani & Ricci 2017; Baloković et al. 2018, 2019; Tanimoto et al. 2019; BuChner et al. 2019, 2021), but either they have a limited geometrical flexibility, or they do not incorporate the full suite of X-ray physics. To our knowledge, there is no radiative transfer code that combines advanced X-ray radiative processes such as bound-electron scattering and interactions with dust grains, with a full-3D treatment of radiative transfer that overcomes geometrical limitations such as the need for 3D transfer media to exhibit certain spatial symmetries, to follow a 2D distribution of spherical clumps, or to be assembled as a superposition of 1D and 2D building-block components.

To this end, we extend the 3D radiative transfer code SKIRT (Camps & Baes 2020) with the X-ray physical processes in gas and dust, to obtain an X-ray radiative transfer code with all benefits of the established SKIRT framework. Beyond the study of circumnuclear media, these X-ray extensions allow for various potential applications of the SKIRT code in X-ray astrophysi-cal research, such as dust mineralogy with XRISM (Corrales et al. 2016; XRISM Science Team 2020), modelling clumpy nova ejecta in supersoft X-ray sources (Ness et al. 2022), studying gas and dust in the interstellar medium (Costantini et al. 2019), constraining AGN contributions in galaxy SED fitting (Yang et al. 2020, 2022), and more.

This work presents the new X-ray functionalities of the SKIRT radiative transfer code, with a particular focus on gas and dust in AGN. The outline for this paper is as follows: in Sect. 2, we introduce the SKIRT radiative transfer code. In Sect. 3, we present the new X-ray processes and their implementation to the SKIRT code. In Sect. 4, we verify our X-ray implementation by comparing SKIRT with some well-established spectral models. Finally, we demonstrate the 3D capabilities of SKIRT in the X-ray regime by calculating images and spectra of clumpy torus models in Sect. 5. In Sect. 6, we summarise our results and discuss their implications.

2 Monte Carlo radiative transfer with SKIRT

Monte Carlo radiative transfer (MCRT) is a direct simulation method to study the effects of absorption, scattering, and re-emission in complex transfer media. In this method, the radiation field is represented as a flow of discrete photon packets, and all relevant physics are emulated by drawing random events from the appropriate probability distributions (Whitney 2011; Noebauer & Sim 2019). The MCRT technique is the most popular numerical method to study radiative transfer in general 3D geometries, partly due to the existence of various acceleration mechanisms that are based on importance sampling (Steinacker et al. 2013). One of these optimisations is the photon peel-off technique to detect flux contributions, which enables efficient simulations in full 3D (Yusef-Zadeh et al. 1984).

SKIRT is a state-of-the-art MCRT code, developed and maintained at Ghent University (Baes et al. 2003, 2011; Camps & Baes 2015, 2020). Historically, the SKIRT code has focussed on dusty astrophysical systems, modelling scattering, absorption, and re-emission by dust. More recently, the SKIRT code was extended to support radiative transfer in electron and gas media, in addition to dust (Camps & Baes 2020). SKIRT offers a large suite of built-in geometries (Baes & Camps 2015), radiation sources, and medium characterisations, in addition to interfaces for post-processing hydrodynamical simulations based on the SPH, AMR, or Voronoi paradigms (Camps et al. 2016; Trčka et al. 2020; Popping et al. 2022). Furthermore, hierarchical octree and k-d tree algorithms are available for efficiently gridding complex transfer geometries (Saftly et al. 2013, 2014). SKIRT has been successfully applied in different fields of astrophysics, with recent applications in galaxies (e.g. Popping et al. 2022; Shen et al. 2022; Vijayan et al. 2022; de Graaff et al. 2022; Camps et al. 2022; Trčka et al. 2022; Jang et al. 2022) and AGN (e.g. Stalevski et al. 2019; Toba et al. 2021; Victoria-Ceballos et al. 2022; Isbell et al. 2022; Stalevski et al. 2023). The SKIRT C++ code is open-source, well-documented¹, and publicly available online², with tutorials for both users and developers.

The latest major code redesign (SKIRT version 9) introduced polarisation by aligned dust grains and line transfer (Camps & Baes 2020; Vandenbroucke et al. 2021; Camps et al. 2021). Furthermore, this SKIRT version supports advanced kinematics, modelling the effect of bulk velocities and velocity dispersions, for both sources and transfer media. in Sect. 3, new X-ray processes are described in the local rest frame of the transfer medium, while SKIRT tracks the flow of photons in the model rest frame (i.e. the model coordinate system). This implies that photon energies are transformed to the moving medium frame before each interaction, and back to the model frame after the interaction: $$E\prime =E\left(1-\frac{k\cdot \upsilon }{c}\right),$$(1)

with υ being the medium velocity, k the photon propagation direction, and E the photon energy in the model frame. In this way, kinematic effects are self-consistently incorporated into the SKIRT radiative transfer results.

SKIRT implements a large number of state-of-the-art acceleration mechanisms to speed up its MCRT calculations. Most notably, the photon peel-off technique (Yusef-Zadeh et al. 1984) allows for recording flux contributions at each scattering event, as opposed to collecting photon packets as they escape the simulated model. In addition to realising a significant speed up, this peel-off technique also avoids the spurious blurring otherwise caused by averaging over a finite solid angle around the observer direction. Other advanced techniques include path length stretching, forced scattering, continuous absorption, composite biasing and explicit absorption (Steinacker et al. 2013; Baes et al. 2016, 2022). Combined, these mechanisms drastically speed up the SKIRT simulations (as compared to other radiative transfer codes), especially in complex 3D geometries. Moreover, the SKIRT code supports a hybrid parallelisation scheme, combining multi-threading and multi-processing, which allows for efficient simulations on systems ranging from laptops to high-performance computing facilities (Verstocken et al. 2017).

Up to this point, the SKIRT code was operational in the millimetre to far-UV wavelength range. in this work, we extend the advanced treatment of radiative transfer into the X-ray range, obtaining an X-ray MCRT code with all benefits of the established SKIRT framework. Furthermore, this allows for radiative transfer simulations that self-consistently cover the X-ray band and the dust-dominated iR-to-UV range, linking X-ray reprocessing to dust modelling in AGN obscuring media.

3 X-ray radiative processes

Focussing on the radiative processes that govern the broad-band spectra of obscured AGN in the 0.1 to 500 keV range, we implement Compton scattering on free electrons, photo-absorption and fluorescence by cold atomic gas, scattering on bound electrons, and extinction by dust. This includes a detailed description of anomalous Rayleigh scattering (Sect. 3.2.3) and extreme-forward scattering by dust (Sect. 3.3.2), which goes beyond the X-ray physics that is currently implemented in AGN obscuration models available for spectral fitting.

3.1 Free-electron media: Compton scattering

Compton scattering describes the inelastic scattering of high-energy photons on free electrons, which forms the relativistic extension of Thomson scattering into the X-ray domain. Comp-ton scattering has a differential scattering cross section given by the Klein & Nishina (1929) formula, which is a function of the incoming photon energy: $$\frac{\text{d}{\sigma }_{\text{KN}}}{\text{dΩ}}\left(\theta ,x\right)=\frac{3{\sigma }_{\text{T}}}{16\pi }\left[C^3\left(\theta ,x\right)+C(\theta ,x)-C^2\left(\theta ,x\right){\sin }^2\theta \right],$$(2)

with x = E/m_ec² the photon energy scaled to the electron rest energy, σ_T ≈ 6.65 × 10⁻²⁵ cm² the total Thomson cross section, and C(θ, x) the Compton factor equal to: $$C\left(\theta ,x\right)=\left(1+x\left(1-\cos \theta \right)\right){.}^{-1}$$(3)

Integration of the differential cross section yields the total Compton scattering cross section as a function of the incoming photon energy: $${\sigma }_{\text{KN}}\left(x\right)=\frac{3{\sigma }_{\text{T}}}{4}\left[\frac{1+x}{{\left(1+2x\right)}^2}+\frac{2}{x^2}+\left(\frac{1}{2x}-\frac{x+1}{x^3}\right)\text{In}\left(2x+1\right)\right].$$(4)

Compared to Thomson scattering, the Compton cross section is reduced at higher energies. However, Compton scattering still forms the dominant source of X-ray extinction above 12 keV; for more details, readers can refer to Sect. 3.2.4. By normalising the Klein-Nishina formula Eq. (2), we obtain the Compton scattering phase function: $$\text{Φ}\left(\theta ,x\right)=\frac{1}{{\sigma }_{\text{KN}}\left(x\right)}\frac{\text{d}{\sigma }_{\text{KN}}}{\text{dΩ}}\left(\theta ,x\right).$$(5)

The energy dependence of this phase function is significant above 10 keV and should not be neglected in detailed radiative transfer simulations. Equation (5) predicts a bias towards forward scattering that increases with increasing photon energy. Even then, Compton scattering is not restricted to forward directions only, and the exact phase function shape should be implemented.

Unlike Thomson scattering, photon energies are changed by Compton scattering on free electrons. Conserving the four-momentum of the photon-electron pair, one has: $$E_{\text{out}}=C\left(\theta ,x\right)E_{\text{in}},$$(6)

with C(θ, x) being the Compton factor as given by Eq. (3). in the electron rest frame, photons can only lose energy upon interacting with a free electron. The relative energy loss is most significant for high-energy photons, and increases with increasing scattering angle. Forward scattering does not change the photon energy, while backscattering maximally reduces the energy by a factor of 1 + 2x. In result, the observed X-ray spectrum at a certain photon energy is formed by photon interactions at higher energies, which must be included in the radiative transfer simulations to obtain accurate reflection spectra, see Sect. 3.2.4. Below 0.25 keV, the energy shift becomes negligible (<0.1%) and Compton scattering converges to elastic Thomson scattering, with a differential cross section: $$\frac{\text{d}{\sigma }_{\text{T}}}{\text{dΩ}}\left(\theta \right)=\frac{3{\sigma }_{\text{T}}}{16\pi }\left[1+{\cos }^2\theta \right],$$(7)

which does not depend on the incoming photon energy.

In SKIRT, the spatial distribution of free electrons is described by an electron number density n_e(r), which defines the transfer geometry of the free-electron medium. We introduce Compton scattering to the SKIRT code using the exact formulae presented in this section. The scattering phase function (Eq. (5)) is an intricate function of the photon energy, which cannot be inverted analytically to generate random scattering angles (Devroye 1986). However, this phase function can be sampled efficiently by using a variation on Khan’s technique, as described by Hua (1997). This algorithm combines the composition and rejection methods to avoid computationally expensive operations, with a rejection rate of about 30% depending on the photon energy.

3.2 Cold-gas media

3.2.1 Photo-absorption

In this work, photo-absorption describes the absorption of X-ray photons in cold atomic gas, where the photon energy is used to liberate a bound electron from a neutral gas atom. Photo-absorption forms the main source of X-ray extinction below 12 keV (see Sect. 3.2.4), but drops in importance at higher energies. This causes Compton-thick AGN to have most of their soft X-ray emission obscured up to some energy, which depends on the total amount of intervening gas.

The photo-absorption cross section σ_Z(E) of a given atom Z is just the sum of its electronic subshell contributions: $${\sigma }_Z\left(E\right)={\displaystyle \sum _{nl}{\sigma }_{Z,nl}}\left(E\right),$$(8)

which displays spectral features reflecting the internal electron structure. Photo-absorption cross sections σ_Z,nl(E) for individual electronic subshells nl are presented by Verner & Yakovlev (1995), who provide analytic fits³ to the cross sections of all ground-state subshells of the first 30 neutral elements (H to Zn), based on Hartree-Dirac-Slater calculations. Updated cross sections for the outer shells of 18 selected elements⁴ are provided by Verner et al. (1996) based on R-matrix calculations, which are valid up to some maximum energy Emax. This update is most important for neutral He, having just one electron K-shell, which dominates the neutral gas extinction in the 0.1 to 0.5 keV range. For neutral hydrogen, both references list the exact same cross section. SKIRT incorporates the Verner & Yakovlev (1995) subshell cross sections for all 30 elements, updated with the available Verner et al. (1996) outer-shell data for energies up to Emax. This implementation is identical to the PHABS model in XSPEC (Arnaud 1996), although SKIRT does include 12 more elements⁵.

Combined with the abundances a_Z for each element, the total photo-absorption cross section per H-atom is calculated as: $${\sigma }_{\text{PA}}\left(E\right)={\displaystyle \sum _Z{a}_Z}\left(1-{\beta }_Z\right)\left({\displaystyle \sum _{nl}{\sigma }_{Z,nl}}\left(E\right)\right),$$(9)

with a_Z being the number density of element Z relative to H I, and 1 – β_Z the fraction of this element found in the gas phase (i.e. not locked in dust grains, see Sect. 3.3). In SKIRT, Anders & Grevesse (1989) solar abundances a_Z are used as the default atomic gas composition. In addition, we allow for custom gas mixes by considering the abundance table as an input parameter. Finally, the spatial distribution of the cold gas medium is characterised by a neutral hydrogen number density n_H(r).

The total photo-absorption cross section for a cold gas medium is shown in Fig. 1, assuming Anders & Grevesse (1989) abundances and no dust depletion. This cross section shows characteristic absorption edges where the individual elements start to contribute to the total gas absorption. Most prominent are the O K-edge at 0.538 keV and the Fe K-edge at 7.124 keV, while the Fe L-edge at 0.724 keV is the only significant L-shell feature.

In general, gas velocity dispersions must be accounted for by convolving all subshell cross sections σ_Z,n1(E) with a Gaussian kernel of appropriate width⁶. In practice however, these cross sections are smooth functions of energy, except for the discontinuity at the absorption edge. Therefore, the effect of a convolution is limited to energies close to the edge energy E_edge, smoothing the sharp step to a sigmoid. This can be incorporated into the subshell cross sections as: $${\sigma }_{Z,nl}^{\prime }\left(E\right)=\frac{{\sigma }_{Z,nl}\left(E_{\text{edge}}+2E_{\text{th}}\right)}{2}\left[1+\text{erf}\left(\frac{E-E_{\text{edge}}}{E_{\text{th}}}\right)\right],$$(10)

for energies between E_edge − 2E_th and E_edge + 2E_th, with E_th being the energy dispersion corresponding to the thermal velocity of the absorbing atom⁷, considering that the convolution of a step function with a Gaussian is given by the error function. This approximation to incorporate velocity dispersions achieves good accuracy and continuity around the edges of all important subshells. However, it is important to note how the neutral atomic material surrounding AGN is expected to have gas temperatures less than 10^3.5 K (see e.g. Wada et al. 2016), which makes the corresponding smoothing unresolvable with current X-ray instrumentation. Nevertheless, the effect of these velocity dispersions can be studied in great detail with SKIRT, which will be especially valuable for studying ionised gas media at higher temperatures in the future.

Figure 1 also shows the contributions of photo-absorption by different electronic shells. In the X-ray range, the dominant absorption process is K-shell photo-absorption, where an inner-shell electron is ejected from the atom. These K-shell electrons have binding energies in the 0.1–10 keV range, where their photo-absorption cross sections peak. Absorption by higher-shell electrons kicks in at lower energies, and quickly becomes more unlikely with increasing energy: σ_Z,nl(E) ∝ E E⁻³. While L-shell absorption is still important for some heavy elements (mainly Fe), the total contribution of M- and N-shells is less than 3.5% over the entire range where photo-absorption is relevant (i.e. below 40 keV, see Sect. 3.2.4). However, including these sub-shells does not slow down our radiative transfer calculations, as the total cross section for a configured gas mix is calculated only once during simulation setup.

Figure 2 shows the relative contributions of each element to the total photo-absorption cross section. Overall, the main contributing elements are just the elements that are most abundant. Below 0.5 keV, the photo-absorption is dominated by helium, and therefore it was important that we incorporated the updated outer-shell cross section of He as provided by Verner et al. (1996). At higher energies, oxygen becomes the most important element, up to 7.1 keV where K-shell absorption by iron kicks in. In spite of its dominant abundance, the contribution of neutral hydrogen is marginal in the X-ray range above 1 keV (<5%).

Currently, SKIRT focusses on photo-absorption by neutral atomic gas, modelling the cold and distant material that causes most of the X-ray extinction in Compton-thick AGN. This same cold gas produces the characteristic X-ray reflection spectrum that is used to study heavily obscured sources, and is possibly located at the distant inner wall of the dusty torus (Sect. 1). In this work, we are not implementing radiative transfer in ionised plasmas, which can be found more closely to the central AGN. However, we verified how the photo-absorption cross sections of mildly ionised material differ only slightly from the neutral cross sections, allowing SKIRT to be used in these types of media. We note that interested users can easily add ionised media to the open-source developing environment of the SKIRT code⁸.

Notes. The associated fluorescent yields and line energies for Fe I are listed as an illustration. We note that Bearden (1967) lists equal line energies for the Fe K_β1 and Fe K_β3 sublines.

Finally, we did not consider neutral elements with Z > 30, which can be neglected because of their much lower abundances.

Fig. 1 Total photo-absorption cross section for a cold gas medium with Anders & Grevesse (1989) abundances and no dust. The K- and L-shells dominate the absorption over the entire X-ray range. The prominent absorption edge structure reflects the atomic composition of the gas.

Fig. 2 Main contributions to the total photo-absorption cross section of a cold gas (see Fig. 1), assuming Anders & Grevesse (1989) abundances and no dust.

3.2.2 Fluorescence

Fluorescence forms one possible way for gas atoms to relax after a photo-absorption event. After an inner-shell electron has been ejected by the absorption of an X-ray photon, a higher-shell electron can fall into this inner-shell vacancy, emitting the energy difference as a fluorescent photon. As fluorescence describes bound-bound transitions in atoms, characteristic line photons are emitted, which produce fluorescent lines in the X-ray spectra of AGN. The most prominent fluorescent line is the Fe K_α line at 6.4 keV, which forms a strong, characteristic feature in the reprocessed X-ray emission of Compton-thick AGN (Mushotzky et al. 1978; Nandra et al. 1989; Pounds et al. 1989). A competing, non-radiative relaxation channel is described by the Auger process, which involves the ejection of one or more electrons, with no effect on the radiation field.

In the X-ray range, we focus on K-shell fluorescence only, which is the radiative relaxation after a K-shell photo-absorption (see Sect. 3.2.1). Table 1 lists the fluorescent transitions that are implemented in SKIRT, with their corresponding line identifiers. K_α-transitions can occur in neutral atoms with Z > 5, while K_β¯transitions are possible for Z > 12. The probability of a fluorescent transition K¡ towards a given K-shell vacancy is expressed by the fluorescent yield $Y_{Z,K_i}$, which is provided by Perkins et al. (1991) for all elements in the photo-absorbing gas⁹. As the L-shell photo-absorption cross sections and fluorescent yields are much lower than the corresponding K-shell values, L-shell fluorescence can be neglected in the X-ray range. Likewise, M- and N-shell fluorescence are negligible relative to the K-shell.

Figure 3 shows the fluorescent yields for all transitions that are incorporated in SKIRT. For each element, the K_αı-transition has the highest probability, producing a fluorescent line that is two times stronger than the corresponding K_α2-line, and 5.5 times stronger than both K_β -lines combined. As the fluorescent yields increase along the atomic sequence, abundant heavy elements like Fe and Ni dominate the fluorescent line emission of AGN. For each transition, the corresponding line energy is taken from Bearden (1967), who provides lab measurements that should be more accurate than the theoretical values as calculated by Perkins et al. (1991). As the K_αı- and K_α2-lines have very similar line energies, they are often observed as a single K_α-line with current CCD instruments. For neutral iron, both sublines are separated by 13 eV (see Table 1), which could be resolved with upcoming microcalorimeter observations (Barret et al. 2018; Tashiro et al. 2020). Analogously, one observes a single K_β-line at slightly higher photon energies.

In principle, fluorescence is a re-emission process that could be performed in the simulation after calculating the total photo-absorption in each gas cell. This implementation would be similar to the dust emission cycle in SKIRT, which needs some iterations due to dust self-absorption. In practice however, the relaxation time after a photo-absorption is short, and the fluorescent re-emission can be considered as being instantaneous (Krause & Oliver 1979). Therefore, we implement fluorescence as a scattering process, where the incoming photon is re-emitted isotropically, with a photon energy updated to the fluorescent line energy. This allows both photo-absorption and fluorescence to be treated during SKIRT’S primary emission cycle (for details, see Camps & Baes 2020).

The effective scattering cross section for fluorescence on an atom Z is just the K-shell photo-absorption cross section σz,1s(E), multiplied by the total fluorescent yield Y_z,ĸ of the K-shell, which is the sum of all four yields shown in Fig. 3. The effective K-shell absorption cross section is then correspondingly updated to (1 – Y_z,_ĸ)σ_z,1s(E), and Eq. (9) should be interpreted as the total extinction cross section due to photo-absorption and fluorescence combined. When fluorescent re-emission occurs in an atom, the scattered photon is assigned a fluorescent line energy sampled from the discrete probability distribution given by the set of yields $\left\{Y_{Z,K_i}\right\}$. Finally, this line energy is Doppler shifted to be consistent with the thermal velocity of the atom (sampled from a Maxwell distribution), in addition to the bulk velocity of the gas, producing self-consistently broadened lines.

Fig. 3 Fluorescent yields for all fluorescent transitions that are implemented in SKIRT, as provided by Perkins et al. (1991). See Table 1 for the transition definitions.

3.2.3 Bound-electron scattering

X-ray scattering in cold atomic gas is mainly caused by the electrons that are bound to the neutral gas atoms. Two different types of bound-electron scattering can be distinguished. Inelastic bound-Compton scattering describes the incoherent, particle-like scattering of X-ray photons on individual bound electrons, while elastic Rayleigh scattering describes the coherent, wave-like scattering on atoms as a whole. Both types of bound-electron scattering are important for modelling the X-ray spectra of obscured AGN, as these interactions make up the reprocessed X-ray continuum, and produce characteristic spectral features such as the Compton hump and the Fe K_α Compton shoulder.

To describe the bound-Compton scattering and Rayleigh scattering processes, we introduce the dimensionless momentum transfer parameter q, following the definition of Hubbell et al. (1975): $$q=\frac{E}{12.4\text{keV}}\sin \left(\frac{\theta }{2}\right),$$(11)

which is a function of the photon energy E and the scattering angle θ. The differential scattering cross section for bound-Compton scattering on atom Z is just the Klein-Nishina formula (Eq. (2)), with a multiplicative correction factor S_z(q): $$\frac{\text{d}{\sigma }_{\text{CS},\text{Z}}}{\text{dΩ}}\left(\theta ,E\right)=\frac{3{\sigma }_{\text{T}}}{16\pi }\left[C^3+C-C^2{\sin }^2\theta \right]S_Z\left(q\right),$$(12)

where C = C(θ, E) is the Compton factor, which does also depend on the scattering angle and the photon energy, see Eq. (3). The correction factor S_z(q) is the incoherent scattering function, which is provided by Hubbell et al. (1975) for all elements in the photo-absorbing gas (see Sect. 3.2.1). This function represents the binding-energy correction to the free-electron description of X-ray scattering (see Sect. 3.1), and is a smooth function of q that converges to Z for large q, and decreases to zero for small q. The total scattering cross section for bound-Compton scattering on element Z can be obtained as: $${\sigma }_{\text{CS},\text{Z}}\left(E\right)=2\pi {\displaystyle {\int }_0^{\pi }\text{d}\theta }\sin \theta \frac{\text{d}{\sigma }_{\text{CS},\text{Z}}}{\text{dΩ}}\left(\theta ,E\right),$$(13)

which needs to be evaluated numerically for each element, based on the tabulated incoherent scattering functions. The total bound-Compton scattering cross section per H-atom is then: $${\sigma }_{\text{CS}}\left(E\right)={\displaystyle \sum _Z{a}_Z}\left(1-{\beta }_Z\right){\sigma }_{\text{CS},\text{Z}}\left(E\right),$$(14)

with a_z being the abundance of element Z relative to H I, andβ_z the correction for atoms locked up in dust grains (see Sect. 3.2.1). The total scattering cross section σ_CS(E) is shown in Fig. 4, assuming Anders & Grevesse (1989) abundances and no dust.

Bound-Compton scattering is most important at high photon energies, where it forms the main source of X-ray extinction, see Sect. 3.2.4. Below 20 keV however, bound-Compton scattering is heavily reduced, as the scattering functions S_z(q) start suppressing the free-electron formula in Eq. (12). The corresponding scattering phase function is: $${\text{Φ}}_{\text{CS}}\left(\theta ,E\right)=\frac{1}{{\sigma }_{\text{CS}}\left(E\right)}{\displaystyle \sum _Z{a}_Z}\left(1-{\beta }_Z\right)\frac{\text{d}{\sigma }_{\text{CS},Z}}{\text{dΩ}}\left(\theta ,E\right),$$(15)

which is illustrated in Fig. 5 for a photon energy of 3.3 keV. Bound-Compton scattering strongly suppresses the forward scattering direction, which was the preferred direction for free-electron scattering, see Sect. 3.1. Finally, bound-Compton scattering is an inelastic process, changing the incoming photon energy as described by Eq. (6).

Analogously, the differential scattering cross section for elastic Rayleigh scattering on element Z is given by the Thomson formula (Eq. (7)), with a multiplicative correction factor $F_Z^2\left(q\right)$ $$\frac{d{\sigma }_{\text{RS},Z}}{\text{dΩ}}\left(\theta ,E\right)=\frac{3{\sigma }_{\text{T}}}{16\pi }\left[1+{\cos }^2\theta \right]F_Z^2\left(q\right),$$(16)

where F_z(q) is the atomic form factor of element Z, which has been calculated by Hubbell et al. (1975) for all elements H to Zn. This form factor represents the quantum mechanical scattering amplitude for X-ray scattering on a cloud of bound electrons, which converges to Z at small q, and decreases to zero at large q. The total cross section for Rayleigh scattering on Z is then: $${\sigma }_{\text{RS},Z}\left(E\right)=2\pi {\displaystyle {\int }_0^{\pi }\text{d}\theta }\sin \theta \frac{d{\sigma }_{\text{RS},Z}}{\text{dΩ}}\left(\theta ,E\right),$$(17)

and from these tables, the total Rayleigh scattering cross section per H-atom can be calculated as: $${\sigma }_{\text{RS}}\left(E\right)={\displaystyle \sum _Z{a}_Z\left(1-{\beta }_Z\right)}{\sigma }_{\text{RS,}Z}\left(E\right),$$(18)

which is shown in Fig. 4. Rayleigh scattering forms the main scattering channel at soft X-ray energies, allowing some soft X-ray photons to escape the system before photo-absorption can occur. Finally, the corresponding scattering phase function can be obtained as: $${\text{Φ}}_{\text{RS}}\left(\theta ,E\right)=\frac{1}{{\sigma }_{\text{RS}}\left(E\right)}{\displaystyle \sum _Z{a}_Z}\left(1-{\beta }_Z\right)\frac{\text{d}{\sigma }_{\text{RS,Z}}}{\text{dΩ}}\left(\theta ,E\right),$$(19)

which is illustrated for 3.3 keV photons in Fig. 5, showing a clear bias towards forward scattering.

A final level of refinement in treating bound-electron scattering is obtained by incorporating the anomalous scattering factors for Rayleigh scattering, which introduce spectral features close to the atomic absorption edges: $$\frac{\text{d}{\sigma }_{\text{RS}\prime ,Z}}{\text{dΩ}}\left(\theta ,E\right)=\frac{3{\sigma }_{\text{T}}}{16\pi }\left[1+{\cos }^2\theta \right]\left[{\left(F_Z\left(q\right)+F{\prime }_Z\left(E\right)\right)}^2+F{″}_Z\left(E\right)\right],$$(20)

with ${F^{\prime }}_Z\left(E\right)$ and ${F^{″}}_Z\left(E\right)$ being the real and imaginary anomalous scattering factors of element Z, which are provided by Cullen et al. (1997) for all 30 elements in the gas. These anomalous scattering factors correct the Rayleigh formula (Eq. (16)) for resonant scattering near the atomic absorption edges, and depend on the photon energy only. The anomalous Rayleigh scattering cross sections of four selected elements are shown in Fig. 6, displaying distinct spectral features at the absorption edge energies. Away from the absorption edges, these cross sections converge to the original description of smooth Rayleigh scattering.

Some X-ray radiative transfer codes are approximating bound-electron scattering by free-electron scattering, which is somewhat easier to implement, see Sect. 3.1. In addition, older codes might have had to consider computational constraints which did not allow for sampling the complex numerical functions that were introduced in this section. This free-electron approximation could be justified above 30 keV, where bound-Compton scattering converges to free-electron scattering, and Rayleigh scattering is negligible, see Sect. 3.2.4. For Anders & Grevesse (1989) solar abundances, the free-electron approximation is equivalent to assuming Compton scattering on 1.21 free electrons per H-atom: $${\displaystyle \sum _Z{a}_{Z}Z=1.21,}$$(21)

which is shown in Figs. 4 and 5 for comparison. However, while this approximation matches the bound-electron scattering cross section where scattering dominates the total gas extinction (i.e. above 30 keV, see Sect. 3.2.4), it underpredicts scattering at lower energies (see Fig. 4), where bound-electron scattering is still important. Furthermore, the free-electron approximation ignores the important subtlety that scattering is elastic at low X-ray energies, and inelastic at high X-ray energies. In addition, this approximation does not incorporate the distinct phase function shapes that govern the X-ray scattering physics (see Fig. 5), and have their effect on the reprocessed X-ray emission.

The computational efficiency of modern radiative transfer codes such as SKIRT combined with the availability of powerful computers now allows for physically more correct simulations that include bound-electron scattering, SKIRT simulations with bound-electron scattering enabled are three times slower than the corresponding free-electron simulations. However, the effect of bound-electron scattering can be observed in the resulting X-ray spectra and should not be neglected, especially in X-ray radiative transfer simulations that cover photon energies below 30 keV. For benchmarking purposes, we provide the options to approximate bound-electron scattering by free-electron scattering, or to fully ignore scattering in cold gas.

Fig. 4 Bound-electron scattering cross section of a cold gas medium, assuming Anders & Grevesse (1989) abundances and no dust. The 1.21 free-electron approximation is shown in grey for comparison (see text).

Fig. 5 Normalised scattering phase functions for bound-electron scattering in cold gas, at a photon energy of 3.3 keV. At this photon energy, bound-Compton scattering and Rayleigh scattering are equally important, see Fig. 4. The free-electron phase function at 3.3 keV is shown in grey for comparison.

Fig. 6 Anomalous Rayleigh scattering cross sections of H, Be, Mg, and Zn. The corresponding smooth Rayleigh cross sections (see Eq. (18)) are shown in black.

Fig. 7 Total gas extinction cross section for Anders & Grevesse (1989) abundances and no dust. The effective photo-absorption (blue), bound-electron scattering (green), and effective fluorescent scattering (yellow) contributions are explicitly shown, to illustrate how these processes shape the X-ray spectrum in radiative transfer simulations.

3.2.4 Total gas extinction

Figure 7 shows the total extinction cross section of cold gas over the 0.1 to 500 keV range. Below 12 keV, the total gas extinction is dominated by photo-absorption (Sect. 3.2.1), which removes most soft X-ray photons from the radiation field. The effective photo-absorption cross section is shown in blue, which is Eq. (9) corrected for fluorescent re-emission as described in Sect. 3.2.2. Fluorescence (yellow) has a limited contribution to the total gas extinction, but shapes the X-ray reflection spectrum by re-emitting photons at fixed line energies (Sect. 3.2.3), producing fluorescent line features.

Above 12 keV, bound-electron scattering (Sect. 3.2.3) starts to dominate the total gas extinction. Figure 7 shows the bound-electron scattering cross section in green, which is the sum of bound-Compton scattering and Rayleigh scattering (see Fig. 4). Bound-Compton scattering forms the main source of cold-gas extinction at high energy, and produces the characteristic X-ray reflection spectrum including the Compton hump and the Fe K_α Compton shoulder. Rayleigh scattering is most important at low X-ray energies, but never dominates the total gas extinction. Yet, it forms the main scattering channel for soft X-ray photons, producing the X-ray reflection component below 4 keV (Fig. 4).

The X-ray reflection continua of AGN circumnuclear media are produced by Compton down-scattering on bound electrons, as described in Sect. 3.2.3. Therefore, the observed spectrum at a given photon energy is formed by photon interactions at higher energies, which have to be included in the radiative transfer simulations. We found that simulated X-ray spectra up to 100 and 150 keV attain full convergence when considering X-ray interactions up to 200 and 400 keV, respectively. In addition, we found that SKIRT simulations can obtain accurate model spectra over the entire spectral range of current hard X-ray observatories such as NuSTAR (Harrison et al. 2013) and Swifi/BAT (Gehrels et al. 2004; Barthelmy et al. 2005; Krimm et al. 2013), and proposed hard X-ray missions such as HEX-P (Madsen et al. 2019).

3.3 Dust media

3.3.1 X-ray dust model

The detailed 3D distribution of dusty gas plays a major role in the IR-to-UV spectral modelling of circumnuclear media, see Sect. 1. Nevertheless, dust has mostly been neglected in X-ray radiative transfer simulations of AGN, even though it could contribute significantly to the total X-ray extinction. Furthermore, the absorption of X-ray photons could alter the dusty medium state, that is heat or destroy the dust grain population, affecting the observations at IR to UV wavelengths. Therefore, we introduce an X-ray dust model to the SKIRT code, extending the advanced treatment of dust radiative transfer into the X-ray range. This will allow for self-consistent model predictions covering both the X-ray band and the IR-to-UV wavelength range, linking X-ray reprocessing to dust modelling in AGN.

As a first X-ray dust model for SKIRT, we implement a dust mix of carbonaceous grains and amorphous silicate grains with a Weingartner & Draine (2001) grain size distribution. This dust model has already been implemented in SKIRT for IR to UV wavelengths (Li & Draine 2001), allowing the dust mix to be used over the entire 10⁻⁴ to 10⁴ µm range. The X-ray extinction cross sections for this model were computed¹⁰ by Draine (2003c) for photon energies between 0.1 keV and 12.4 keV, and represent the absorption and scattering probabilities of individual dust grains that are representative for the entire dust mix. Absorption by dust grains is similar to photo-absorption by the neutral gas atoms that constitute the dust mix composition (Sect. 3.2.1), but accounts for shielding at soft X-ray energies and interference in the solid grains. Scattering on dust grains is an elastic process, which was modelled using Mie theory or anomalous diffraction theory depending on the photon energy (see Draine 2003c, for more details). These cross sections are shown in Fig. 8, and contain detailed spectral features near the atomic absorption edges, which form an interesting diagnostic on the dust composition, crystallinity, shape, and grain size, and have been observed both in astrophysical systems and lab measurements (Corrales et al. 2016; Costantini et al. 2019; Zeegers et al. 2019; Rogantini et al. 2020; Psaradaki et al. 2020). More advanced dust models with configurable size distributions can easily be introduced to SKIRT, owing to the modular design of the code.

Finally, a dust mass density ρ^d(r) characterises the spatial distribution of the dusty medium. Often however, it is more meaningful to connect the distribution of dust to the distribution of the gas in which the dust resides, by declaring a dust mass per H-atom m_d, so that: $${\rho }_{\text{d}}\left(r\right)=m_{\text{d}}{n}_{\text{H}}\left(r\right).$$(22)

For a realistic dust-to-metal mass ratio of 0.3 (Camps et al. 2016) and Anders & Grevesse (1989) solar abundances, this dust mass m_d would be 1.38 × 10⁻²⁶ g per H-atom.

In turn, the atomic gas abundances should be updated to account for the depletion of elements that are locked up in dust grains, see Eq. (9). The Weingartner & Draine (2001) dust mixture implemented in SKIRT consists of 84% graphite and 16% olivine (i.e. the relative number of structural units C and MgFeSiO₄, corresponding to dust mass fractions of 27% and 73%, respectively). Defining the graphite abundance a_gr as the number of structural units per H-atom, one has: $$a_{\text{gr}}=\frac{0.84m_{\text{d}}}{0.84{\mu }_{\text{gr}}+0.16{\mu }_{\text{ol}}},$$(23)

with µ_gr = 12 amu the mass per structural unit graphite, and µ₀₁ = 172 amu the mass per structural unit olivine. The carbon depletion factor is then calculated as β_C = a_gr/a_C. Equivalently, the gas abundances of Mg, Fe, and Si are reduced by the abundance of olivine: $$a_{\text{ol}}=\frac{0.16m_{\text{d}}}{0.84{\mu }_{\text{gr}}+0.16{\mu }_{\text{ol}}},$$(24)

while the oxygen depletion is calculated as β_O = 4a_oı/ao. For solar gas abundances, md should be limited to 1.38 × 10⁻²⁶ g per H-atom, when all available Si would be locked up in dust grains. Phenomenological dust masses exceeding this value could then be explained by the existence of elongated or porous dust grains, which provide more extinction than an equivalent dust mass in compact spherical grains, as implemented in our dust model (Draine 2003a).

The scattering phase function for the X-ray dust model in SKIRT is approximated by the Henyey & Greenstein (1941, HG) function: $${\text{Φ}}_{\text{HG}}\left(\theta ,E\right)=\frac{1}{4\pi }\frac{1-g^2\left(E\right)}{{\left(1+g^2\left(E\right)-2g\left(E\right)\cos \theta \right)}^{3/2}},$$(25)

which has also been used to model dust scattering in the IR-to-UV wavelength range (Draine 2003b). The Φ_HG function shape is set by an energy-dependent asymmetry parameter g(E), which can be interpreted as the expected value for the cosine of the scattering angle: −1 ≤ g = 〈cosθ〉 ≤ 1. In the X-ray range, we adopt a value for g in function of the photon energy E: $$g\left(E\right)=1-\frac{1.23\times {10}^{-3}\text{kev}}{E},$$(26)

which was obtained by fitting a HG function to the dust scattering phase function Φ_D03 proposed by Draine (2003c), for each photon energy. The resulting functional form g(E) then matches the tabulated g values adopted by SKIRT in the millimetre to far-UV wavelength range at 0.1 keV, and causes the Φ_HG and Φ_D03 phase functions to coincide at θ = 0 for all photon energies.

X-ray scattering by dust can be observed directly in diffuse scattering halos around galactic X-ray point sources, which form a promising diagnostic on the dust composition and chemistry (Draine & Tan 2003; Corrales & Paerels 2015; Corrales et al. 2017; Sguera et al. 2020; Lamer et al. 2021; Landstorfer et al. 2022). These dust scattering halos have all been modelled with simple ID dust screens (see e.g. Costantini & Corrales 2022). With SKIRT, more complex dust geometries can now be explored, and tested against X-ray observations.

Fig. 8 Dust extinction cross sections in the X-ray range, as implemented in SKIRT. The atomic absorption edges of the elements that make up the dust grains are indicated in grey.

3.3.2 Extreme forward scattering by dust

As noted in Sect. 3.3.1, the current X-ray dust model in SKIRT describes scattering on dust grains using the HG phase function Φ_HG (Eq. (25)), with asymmetry parameter values g(E) given by Eq. (26). From the latter equation, it is clear that g closely approaches unity in the X-ray regime, causing the corresponding HG function to model extreme forward scattering, with 〈cos θ〉 = 0.9988 at 1 keV, and 〈cos 9〉 = 0.9999 at 10 keV. Evaluating the HG function or sampling scattering angles from this function becomes numerically unstable for large g > 1–10⁻⁶. However, these values of g are only reached for energies E > 1 MeV, well above the considered simulation domain. Nevertheless, the SKIRT implementation clips any g value to the stated numerical limit, to avoid numerical instabilities with future dust mixes that could reach this limit at a lower photon energy.

More problematically, this extreme forward scattering causes issues with the scattering peel-off technique, which forms a core element of SKIRT’S radiative transfer method, as described in Baes et al. (2011). Specifically, the bias factor w_0bS for peel-off photon packets sent to an observer becomes extremely large when the scattering angle towards that observer is small, as the sharply peaked HG function yields large values for small angles; for more details, readers can refer to Eq. (30) of Steinacker et al. (2013), for example. As a result, some peel-off photon packets carry an extreme, artificially boosted amount of energy, which causes unacceptable noise levels in the fluxes measured by SKIRT instruments’ (which are hereafter defined as apertures recording photon packets leaving the system in a preset observer direction).

The bias factor for HG scattering relative to the bias factor for isotropic scattering is: $$b\left(\theta ,E\right)=\frac{w_{\text{obs,HG}}\left(\theta ,E\right)}{w_{\text{obs,iso}}}=\frac{{\text{Φ}}_{\text{HG}}\left(\theta ,E\right)}{{\text{Φ}}_{\text{iso}}}=4\pi {\text{Φ}}_{\text{HG}}\left(\theta ,E\right),$$(27)

which equals a renormalised form of the HG function. For a given photon energy, the largest possible value of b corresponds to peel-off photon packets emitted to an observer in the original photon propagation direction, so that θ_obs = 0. Experiments show that b should be kept under ≈10³, to achieve acceptable noise levels in SKIRT simulations with a manageable number of photon packets. However, for photon energies of 12.4 keV, we find b(θ = 0) ≈ 2.0 × 10⁸, which is many orders of magnitude higher. Even for 0.1 keV photons, we find b(θ = 0) ≈ 1.3 × 10⁴. Thus, this peel-off problem exists for dust scattering over the entire X-ray range under consideration.

To address this situation, we average the phase function part of the bias factor ω_obs over a small patch of sky around the observer direction k_obs, as opposed to evaluating the phase function Φ_HG at θ_obs exactly. More precisely, we replace b by its average value 〈b〉 over a small solid angle Ω_obs around k_obs¹¹: $$\langle b\rangle \left({\theta }_{\text{obs}},E,\delta \right)=\frac{{\displaystyle {\int }_{\cos \left({\theta }_{\text{obs}}+\delta \right)}^{\cos \left({\theta }_{\text{obs}}-\delta \right)}4\pi {\text{Φ}}_{\text{HG}}\left(\cos \theta ,E\right)\text{\hspace{0.17em}d\hspace{0.17em}cos\hspace{0.17em}}\theta }}{\cos \left({\theta }_{\text{obs}}-\delta \right)-\cos \left({\theta }_{\text{obs}}+\delta \right)},$$(28)

where θ_obs is the angle between the original photon direction and the observing instrument, and 6 can be interpreted as the half-opening angle of the instrument in the polar direction. Averaging in the azimuthal direction has no net effect, as the HG phase function does not depend on ϕ.

In the polar direction, we adopt a parameter value of δ = 4°. Indeed, this value for 6 dampens the relative bias factor to 〈b〉 ≲ 820 for all E, comfortably below the empirical limit of ≈10³ that was stated above. Smaller δ values would still yield bias factors at or above the empirical limit, causing unacceptable noise levels in the recorded fluxes. On the other hand, larger 6 values would unnecessarily increase the artificial spatial blur in the results (see next paragraph). For the selected value of δ = 4^°, the averaged relative bias factor 〈b〉 is shown for different photon energies in Fig. 9, which is strongly dampened for small scattering angles, and converges to its original value for large scattering angles. The expression for averaging the bias factor (Eq. (28)) can be evaluated analytically at each scattering peel-off. Therefore, this implementation of extreme forward scattering is equally efficient as the original treatment without averaging.

While decreasing the noise levels as desired, this averaging process also introduces a certain amount of unwanted spatial blur to the recorded fluxes. Indeed, photons moving in a direction close to the line of sight of some instrument are smeared out over a small solid angle, when scattering on a dust grain towards that instrument. However, this blurring is not really an issue when comparing simulation results to spatially resolved observations of dust scattering, as only the scattered flux is affected, mainly in the most central regions¹² of extended sources that are probably dominated by the direct source emission anyway. Furthermore, these central regions do also feature some spatial blur in actual observations, as described by the point spread function. In any case, SKIRT properly recovers the extended halo emission away from the source centre, which can be used to model halo profiles in dust studies. For X-ray sources that cannot be spatially resolved (e.g. AGN), the effect of spatial averaging is negligible, as the observed fluxes already entail a spatial integration over the source. Finally, we recall that the averaging procedure is only invoked at the last scattering event before a photon detection (i.e. at the scattering peel-off), and that all other photons performing a random walk inside the dusty model geometry are treated exactly, based on the unmodified, sharply peaked HG function.

Fig. 9 Dampened relative bias factors 〈b〉 for dealing with extreme forward dust scattering. For δ = 4°, we maintain 〈b〉 ≲ 820.

4 Benchmark tests

To verify our X-ray implementation, we compare the SKIRT simulation output¹³ with some well-established spectral models that deal with X-ray obscuration and reflection in the context of AGN. In Sect. 4.1, we reproduce the setup of some 1D XSPEC models for free-electron, cold-gas, and dust extinction. In Sect. 4.2, we benchmark our code against some popular smooth torus models that are frequently used to model observational X-ray data. In all simulations, the primary X-ray source has a power law shape with a photon index of Γ = 1.8, an exponential cut-off at E_cut = 200 keV, and an integrated 2–10 keV luminosity of L_2–10 = 2.8 × 10⁴² erg s⁻¹, which is representative for the nearby AGN in the Circinus galaxy (see e.g. Wada et al. 2016).

4.1 One-dimensional XSPEC models

4.1.1 PHABS and CABS

PHABS and CABS are two 1D XSPEC models, dealing with cold-gas and free-electron extinction (Arnaud 1996). PHABS models photo-absorption (Sect. 3.2.1) by a 1D slab of neutral gas atoms, evaluating the 1D transmission function: $$T\left(E\right)={\text{e}}^{-N_{\text{H}}{\sigma }_{\text{PA}}\left(E\right)},$$(29)

with N_H being the hydrogen column density of the slab, and σ_PA(E) the photo-absorption cross section per H-atom (Eq. (9)). SKIRT implements the exact same photo-absorption cross sections as those set by the XSPEC command XSECT VERN, allowing for a direct comparison. In addition, we can choose the same gas abundances as those set by the XSPEC command ABUND ANGR, which is Anders & Grevesse (1989) abundances for all 18 elements implemented in PHABS. The PHABS model deals with photo-absorption only, ignoring fluorescence (Sect. 3.2.2) and bound-electron scattering (Sect. 3.2.3) in the gas.

CABS models the extinction of a 1D slab of free electrons (Sect. 3.1), ignoring the actual Compton-scattered photons. This simplified treatment is valid in optically thin media, and comes down to evaluating the ID transmission function: $$T\left(E\right)={\text{e}}^{-N_{\text{e}}{\sigma }_{\text{KN}}\left(E\right)},$$(30)

with N_e being the electron column density of the slab, and σ_KN(E)the Compton scattering cross section (Eq. (4)). The electron column density of the CABS model is often related to the hydrogen column density of an other PHABS model as N_e = 1.21 × N_H, forming a crude approximation to bound-electron scattering in that cold gas medium, see Sect. 3.2.3. In this subsection, we combine the PHABS and CABS models to represent this approximation.

We mimic the XSPEC model setup in SKIRT with a 3D slab of cold gas and free electrons, so that an X-ray point source is observed through the slab. We disable bound-electron scattering, which is a configurable option of the cold-gas mix in SKIRT (see Sect. 3.2.3). We run SKIRT simulations for different values of N_H, each time linking the electron column density to the hydrogen column density as N_e = 1.21 × N_H. As the PHABS and CABS models ignore both fluorescent line emission and scattered photons, we focus on the direct observed flux only (i.e. transmitted photons which did not interact within the slab). It is possible to retrieve this individual flux component using the ‘smart’ photon instruments in SKIRT, which can access in-simulation information such as the interaction history of detected photons, which is not available in actual observations.

The SKIRT simulation output is shown in Fig. 10, in excellent agreement with the corresponding XSPEC results. And even though the PHABS and CABS models do not incorporate any non-trivial radiative transfer effects, they assure a correct implementation of the basic radiative processes in free-electron and cold-gas media. Equivalently, we ran SKIRT simulations with cold gas only, this time enabling the free-electron approximation for bound-electron scattering in the gas (see Sect. 3.2.3). In these simulations, we obtained the exact same results as the ones shown in Fig. 10, in agreement with the underlying assumption of scattering on 1.21 free electrons per H-atom.

Fig. 10 First SKIRT benchmark (Sect. 4.1.1), comparing three SKIRT simulations with the PHABS and CABS models for cold-gas and free-electron extinction. The hydrogen column density NH is indicated on the figure, while the electron column density isNe = 1.21 × NH (see text).

4.1.2 PEXMON 

PEXMON (Nandra et al. 2007) is the standard XSPEC model for dealing with reflection by neutral material, which combines the PEXRAV (Magdziarz & Zdziarski 1995) continuum reflection model with self-consistently generated Fe K_α, Fe K_ß, and Ni K_α fluorescent lines, plus a Fe K_α Compton shoulder. For a PEXMON model parameter REL_REFL = −1, this model assumes a semi-infinite slab of optically thick atomic gas, illuminated by an X-ray point source. The reflected X-ray flux is then observed from the same side of the slab, for different orientations with respect to the reflecting plane. Inside the gas, all H and He atoms are assumed to be fully ionised, and Compton scattering (Sect. 3.1) on these free electrons produces the observed reflection continuum. All other elements are assumed to be neutral, causing photo-absorption similar to the PHABS model (Sect. 4.1.1). In addition, three fluorescent lines are incorporated into the model (Fe K_α, Fe K_ß, and Ni K_α), based on calculations by George & Fabian (1991). Bound-electron scattering is ignored for all elements in the gas.

We mimic the PEXMON model setup in SKIRT using a 3D slab of cold gas and free electrons, with an X-ray point source close to the slab, and an observer at the same side. We set Anders & Grevesse (1989) abundances for all elements in the photo-absorbing gas, except for H and He which are set to zero. Furthermore, we disable any kind of scattering by the gas. We link the density of the free-electron medium to the cold gas as n_e = 1.20 × n_H, which corresponds to a free-electron population formed by ionising all H and He atoms in a Anders & Grevesse (1989) gas mix. We set the hydrogen column density¹⁴ through the slab to an arbitrary high value of N_H = 10²⁷ cm⁻², and the corresponding electron column density to N_e = 1.20 × 10²⁷ cm⁻². These high column densities result in an optical depth through the slab of τ > 500 over the entire simulation domain, which makes the slab virtually impenetrable, modelling an optically thick semi-infinite slab of cold material.

The total observed flux for this model consists of a direct source component, plus a contribution of scattered photons (i.e. photons which interacted at least once with the transfer medium). In this subsection, we focus on the scattered component only, as only scattered photons can reveal the spectral effect induced by radiative transfer processes. Spectra of scattered photons can be retrieved as separate components with both PEXMON and SKIRT, which can be compared to benchmark the SKIRT implementation. We note how fluorescent line photons are considered to be part of the scattered flux component, adding lines to the X-ray reflection spectrum.

Figure 11 shows the SKIRT simulation results for the PEXMON model setup, for three different observing angles i relative to the reflection-plane normal. The SKIRT output spectra show many more fluorescent lines, which is expected as SKIRT implements four line transitions for each element in the photo-absorbing gas (Sect. 3.2.2). The reflected X-ray continuum is found to be in excellent agreement with the PEXMON model over the entire simulation domain, except for the 12–60 keV range, where PEXMON overestimates the scattered flux by up to 10%. However, this behaviour of the underlying PEXRAV reflection model was identified by Magdziarz & Zdziarski (1995), and is therefore expected.

Fig. 11 Second SKIRT benchmark (Sect. 4.1.2), comparing three SKIRT simulations with the PEXMON model for reflection by neutral material. The inclination relative to the reflecting plane is indicated on the figure. The 10% discrepancy in the 12–60 keV range was expected (see text).

Fig. 12 Third SKIRT benchmark (Sect. 4.1.3), comparing SKIRT with the ISMDUST model for dust extinction. The dust mass column density Σd is indicated on the figure. Both models are in good agreement above the C K-edge at 0.280 keV (see text).

4.1.3 ISMDUST 

ISMDUST is a 1D XSPEC model for dust extinction, modelling a 1D slab of graphite and olivine dust grains (Corrales et al. 2016). This dust model is based on the same optical dust properties as those incorporated in SKIRT, but it assumes a different, power law grain size distribution. The ratio of graphite over olivine dust can be varied by tuning the corresponding mass column densities Σ_gr and Σ_ol. In particular, the mass ratio of the dust model incorporated in SKIRT (Sect. 3.3.1) is recovered for: $$\begin{array}{l}{\sum }_{\text{gr}}=0.27\times {\sum }_{\text{d}}\hfill \\ {\sum }_{\text{ol}}=0.73\times {\sum }_{\text{d}},\hfill \end{array}$$(31)

with Σ_d being the total mass column density of the dust.

The ISMDUST model setup is mimicked in SKIRT with a 3D slab of dust between the source and the observer, similar to the geometrical setup of the PHABS and CABS model in Sect. 4.1.1. The corresponding SKIRT simulation results are shown in Fig. 12, for three different dust mass column densities Σ_d. When comparing the direct photon fluxes, the ISMDUST model spectra are nicely recovered for photon energies above 0.29 keV, showing detailed spectral features near the atomic absorption edges of C (0.29 keV), O (0.54 keV), Fe (0.72 keV), Mg (1.31 keV), and Si (1.85 keV). A minor (<8%) offset can be observed in the continuum flux between 0.6 and 3 keV, which could be attributed to the difference in grain size distribution.

Below 0.29 keV (i.e. the C K-edge), the SKIRT and ISMDUST spectra appear to diverge. However, we note that the ISMDUST model was introduced to XSPEC for fitting observational X-ray data from the Chandra and XMM-Newton observatories, and was not designed to be applied below 0.3 keV. For example, the grain size distribution does not include any small grains with sizes < 5 nm, which would increase the extinction at the lowest X-ray energies. The SKIRT dust model on the other hand, is well-defined below 0.3 keV, and extends down to cover the entire X-ray to millimetre wavelength range.

4.2 Radiative transfer torus models

4.2.1 MYTORUS

MYTORUS (Murphy & Yaqoob 2009) is an X-ray spectral model for absorption and reflection by a toroidal reprocessor of cold atomic gas that pioneered the X-ray torus modelling field by providing the first XSPEC model for observational data fitting¹⁵. MYTORUS models a ring torus (i.e. a doughnut) geometry of gas, with a fixed opening angle¹⁶ of 60° and a uniform density. photo-absorption is implemented based on Verner & Yakovlev (1995) and Verner et al. (1996) cross sections and Anders & Grevesse (1989) abundances (i.e. similar to SKIRT, see Sect. 3.2.1), while fluorescence is limited to the iron K_α and K_ß lines, and bound-electron scattering is approximated as free-electron scattering (see Sect. 3.2.3). The primary X-ray source has a power law spectrum similar to the spectrum defined in Sect. 4, but extends up to 500 keV without exponential cut-off.

We reproduce the MYTORUS model in SKIRT with a smooth ring torus of cold gas, centred around an X-ray point source. We ignore bound-electron scattering, and assume Anders & Grevesse (1989) abundances. We run SKIRT simulations for two viewing angles, so that one line of sight sees the central source unobscured (i = 45^°), while the other one intersects the torus (i = 75^°). We vary the equatorial hydrogen column density of the torus between 10²² and 10²⁵ cm⁻², and use the MYTORUS power law input spectrum as described above.

Figure 13 shows the SKIRT simulation results together with the MYTORUS reference spectra. For this comparison, we focus on the reprocessed emission only, which is more informative than the total observed spectra that include direct emission (which can be obtained analytically¹⁷, without any nonlinear radiative transfer effects). We find an excellent agreement between both codes over the entire simulated energy range, for both sightlines and all torus densities. The common Fe K_α and K_ß lines are matched perfectly, as well as the Fe K_α Compton shoulder, while SKIRT implements a large number of additional fluorescent lines. The SKIRT results also reproduce the ‘knee’ feature in the MYTORUS spectra at 169 keV, which is caused by the sudden lack of input photons beyond 500 keV to seed the Compton down-scattered continuum above 169 keV.

For high column densities, one can observe a small offset between the SKIRT and MYTORUS results at low photon energies. However, one should note that column densities of 10²⁵ cm⁻² correspond to optical depths as high as τ > 100 and τ > 2000 for photon energies below 3 and 1 keV, respectively (see Fig. 7). Therefore, one cannot expect realistic soft X-ray model spectra for these heavily obscured sightlines (see e.g. Camps & Baes 2018). Finally, some differences should be expected as the SKIRT spectra are calculated for an observer at the exact inclination of i = 45° and i = 75°, while the MYTORUS model spectra where obtained by averaging the observed radiation over large angular bins centred on 12 representative inclination values, which then need to be interpolated by XSPEC.

Fig. 13 Fourth SKIRT benchmark (Sect. 4.2.1), comparing SKIRT simulations of a torus model against the corresponding MYTORUS model spectra. The opening angle of the ring torus is 60°. The equatorial hydrogen column density NH of the torus is indicated on the figure, for each simulation. Left: scattered photon flux, for a line of sight not intersecting the obscuring torus (i = 45°). Right: scattered photon flux, looking through the obscuring torus (i = 75°). We find an excellent agreement between both codes, for all column densities and sightlines.

4.2.2 REFLEX

REFLEX (Paltani & Ricci 2017) is a well-established ray-tracing code, focussing on X-ray radiative transfer in the circumnuclear media of AGN. The code implements a complete set of X-ray physics in cold-gas media, including bound-electron scattering and a large collection of fluorescent line transitions. In particular, the REFLEX code implements the same interaction cross sections as those that are incorporated in SKIRT, forming an ideal reference for benchmarking the new X-ray processes.

RXTORUS is a smooth torus model that was calculated with REFLEX, and is publicly available for spectral fitting¹⁸. The RXTORUS model represents a ring torus of cold gas, with a uniform density and a variable covering factor, which models scattering on bound electrons (Sect. 3.2.3), in addition to photo-absorption (Sect. 3.2.1) and fluorescence (Sect. 3.2.2). For a torus opening angle of 60° (corresponding to an axis ratio of r/R = 0.5), the RXTORUS geometry is identical to the MYTORUS model geometry (Sect. 4.2.1). RXTORUS has been tested against preceding torus models such as MYTORUS (Murphy & Yaqoob 2009), BNTORUS (Brightman & Nandra 2011), and CTORUS (Liu & Li 2014, 2015), and was applied to the observational X-ray data of NGC 424 (Paltani & Ricci 2017).

We reproduce the RXTORUS model in SKIRT with a uniform ring torus of cold gas, centred around an X-ray point source. Consistent with the RXTORUS model, we enable bound-electron scattering, and assume Anders & Grevesse (1989) abundances. We fix the torus opening angle to 60°, and run SKIRT simulations for an unobscured (i = 45°) and an obscured (i = 75°) sightline. The equatorial hydrogen column density of the torus is varied between 10²² and 10²⁵ cm⁻², which are the N_H-limits of the RXTORUS model.

The SKIRT simulation results are shown in Fig. 14, together with the REFLEX reference spectra. Both models are found to be in excellent agreement over the entire simulation domain, for both sightlines and all column densities between 10²² and 10²⁵ cm⁻². The main observed difference is the high level of photon counting noise in the RXTORUS spectra, especially when the reflected flux is low, for instance when the photo-absorption is high (E < 1 keV) or when there is little reflecting material (N_H = 10²² cm⁻²). The Poisson noise in the corresponding SKIRT results is significantly lower, partly due to the MCRT optimisation mechanisms that were described in Sect. 2. Furthermore, the simulated SKIRT spectra have a higher spectral resolution, producing more narrow fluorescent lines and smoother reflection continua. However, both the fluorescent line fluxes and the continuum levels are found to be consistent in both radiative transfer calculations.

Figure 15 shows a zoom into the important 6.0 to 7.5 keV spectral region, for the same ring torus model with N_H = 10²⁴ cm⁻² and i = 45°. This range contains the three most prominent fluorescent lines (Fe K_α, Fe K_ß, and Ni K_α), plus the Fe K_α Compton shoulder, which forms a powerful probe on the circumnuclear material, being sensitive to both the medium geometry and the physical state of the scattering electrons (bound or free). Despite the limited spectral resolution of the REFLEX simulations, we distinguish clear Compton shoulders in both simulations, with similar strengths and consistent spectral shapes. The SKIRT spectrum reveals some additional substructure in the Compton shoulder shape, which could be observed with upcoming microcalorimeter observatories. Due to the difference in spectral resolution between SKIRT and REFLEX, it is difficult to compare the fluorescent line strengths by eye. Therefore, we measure the integrated line fluxes for Fe K_α, Fe K_ß, and Ni K_α, which we have listed in Table 2. We find an excellent agreement between the line fluxes in both codes, with relative differences in the iron line fluxes as low as 0.9% (Fe K_α) and 0.3% (Fe K_ß), while the much fainter Ni K_α line deviates only 4.5% in flux. Finally, we note how the superior spectral resolution of the SKIRT spectra reveals the K_αı and K_α2 sublines, as shown in Fig. 15.

Notes. The line fluxes in both radiative transfer calculations are found to be in excellent agreement.

Fig. 14 Fifth SKIRT benchmark (Sect. 4.2.2), comparing SKIRT simulations of a torus model against the corresponding RXTORUS model spectra as calculated with REFLEX. The opening angle of the ring torus is 60°. The equatorial hydrogen column density NH of the torus is indicated on the figure, for each simulation. Left: scattered photon flux, for a line of sight not intersecting the obscuring torus (i = 45°). Right: scattered photon flux, looking through the obscuring torus (i = 75°). We find an excellent agreement between both codes, for all column densities and sightlines.

Fig. 15 Zoom on the 6.0 to 7.5 keV reflection spectrum of a ring torus model with NH = 1024 cm−2 (same model as Fig. 14, left). Table 2 lists the integrated line fluxes for the three most prominent fluorescent lines. We note how SKIRT can spectrally resolve the Kα1 and Kα2 sublines.

4.2.3 BORUS

BORUS (Baloković et al. 2018, 2019) is a recent suite of X-ray spectral models for fitting observational data of obscured AGN. Their most popular BORUS02 model implements a uniform-density sphere with conical cutouts at the poles, forming a wedge torus (i.e. a flared disk geometry). This torus geometry is identical to the BNTORUS (Brightman & Nandra 2011) model, and is publicly available as an XSPEC table model¹⁹. BORUS02 has similar model parameters as the RXTORUS model described in Sect. 4.2.2, with additional control over the exponential cutoff energy and the Fe abundance. However, BORUS02 is more restricted in terms of X-ray physics, as it ignores bound-electron scattering, and implements the NIST/xcoM²⁰ photo-absorption cross sections that are limited to E ≥ 1 keV. The BORUS 02 torus geometry forms the most advanced transfer geometry within BORUS, with one more spectral model representing a uniform-density sphere.

We reproduce the BORUS02²¹ model in SKIRT with a uniform wedge torus of cold gas, assuming Anders & Grevesse (1989) gas abundances and free-electron scattering. We fix the torus opening angle to 60°, corresponding to a covering factor of 50%. We run SKIRT simulations for unobscured (cos i = 0.75) and obscured (cos i = 0.25) sightlines, and vary the equatorial hydrogen column density of the torus between 10²² and 10²⁴ cm⁻². These parameter combinations are grid points of the BORUS02 table model, for which the BORUS radiative transfer code was previously run. in this way, we assure that our benchmark results are unaffected by interpolation effects between gridded BORUS02 realisations. Figure 16 shows the SKIRT simulation results and the BORUS reference spectra. The BORUS02 model provides scattered photon fluxes without the direct emission component, which allows for a direct comparison of the X-ray reflection spectra as calculated with SKIRT and BORUS.

We find significant discrepancies between the results obtained with SKIRT and BORUS that are not really understood (see Fig. 16). First, we observe a mismatch between the SKIRT and BORUS reflection spectra for E > 30 keV, which is present for all torus densities and both viewing angles. in this energy range, the reflected X-ray emission is governed by Compton scattering, with photo-absorption being negligible (Fig. 7). Yet, we find the Compton-reflected SKIRT continua to be 20 to 70% lower than the corresponding BORUS results in the 30 to 300 keV range. These high energy differences where traced back to a BORUS02 error in the convolution of the BORUS Green functions with the input spectrum, which will be corrected in an upcoming BORUS02 table release (Balokovic, priv. comm.).

Furthermore, we find that the SKIRT fluxes at lower energies (E < 10 keV) are also disagreeing with the BORUS results, although for unobscured sightlines only (cos i = 0.75, as shown in the left panels of Fig. 16). In this energy range, the reflected X-ray continuum is produced by Compton scattering, but is then heavily photo-absorbed by the intervening gas. Remarkably, we do find a good agreement between both codes for obscured sight-lines in the same energy range (Fig. 16, right). Finally, we obtain an excellent agreement over the 6.0 to 7.5 keV spectral range for obscured sightlines, reproducing the characteristic spectral shape of heavily obscured AGN with both codes (Fig. 16, right).

Our simulated SKIRT spectra have a much higher spectral resolution than the publicly available BORUS02 tables, as needed to model the highest resolution CCD-based spectra and upcoming microcalorimeter observations with XRISM. Even with this larger number of spectral bins, however, our SKIRT results still exhibit less numerical noise.

Fig. 16 Sixth SKIRT benchmark (Sect. 4.2.3), comparing SKIRT simulations of a torus model against the corresponding BORUS02 model spectra as calculated with BORUS. The opening angle of the wedge torus is 60°. The equatorial hydrogen column density NH of the torus is indicated on the figure, for each simulation. Left: scattered photon flux, for a line of sight not intersecting the obscuring torus (cos i = 0.75). Right: scattered photon flux, looking through the obscuring torus (cos i = 0.25). Above 30 keV, we find a significant mismatch between both codes for all sightlines. For unobscured sightlines (cos i = 0.75), the simulation results also diverge for energies below 10 keV.

4.2.4 XARS 

Most recently, the Python-based XARS code (2019) introduced several new medium geometries to X-ray radiative transfer simulations. XARS is an MCRT code designed to study obscuration variability in AGN, which applies the same simulation method as BNTORUS (Brightman & Nandra 2011). The XARS code can handle non-axisymmetric transfer media such as warped disks, clumpy tori, and arbitrary constellations of spherical clumps, and was used to compare hydrodynamical models against observational data of the Circinus galaxy (Buchner et al. 2021). The XARS code is open-source and publicly available online²², and can easily be extended with new transfer geometries.

On the other hand, this Monte Carlo code does not implement some of the well-known biasing techniques that make MCRT simulations efficient in 3D. For example, XARS allows photons to escape the model in all directions, and then averages their contributions over direction bins to obtain the observed flux. Such a detection scheme is inefficient compared to the photon peel-off method (see Sect. 2), especially for non-axisymmetric transfer media. Compared to the REFLEX code (Sect. 4.2.2), XARS does not implement bound-electron scattering, while fluorescence is limited to 11 line transitions, based on fluorescent yields by Bambynek et al. (1972).

We compare SKIRT against the XARS code, focussing on the WEDGE-CUTOFF torus model, which represents a uniform torus with a variable opening angle similar to the BORUS 02 model (see Sect. 4.2.3). We reproduce the model in SKIRT with a uniform wedge torus of cold gas with an opening angle of 60°, identical to the geometrical setup in Sect. 4.2.3. We assume free-electron scattering and Anders & Grevesse (1989) gas abundances, and vary the equatorial hydrogen column density between 10²¹ and 10²⁴ cm⁻². The SKIRT simulation results are shown in Fig. 17, together with the XARS reference spectra for two different viewing angles (45° and 75°). As the WEDGE-CUTOFF model does not provide the scattered flux as a separate model component, we compare the total observed fluxes as calculated with both codes, which causes the reflection spectrum to be hidden by the strong direct component for N_H < 10²³ cm⁻².

We find a good agreement between the SKIRT and XARS results, for both sightlines and all torus densities. The main observed difference is the slightly higher degree of photo-absorption in the XARS spectra at soft X-ray energies, as shown in the right panels of Fig. 17. This offset can be traced back to the photo-absorption cross sections implemented in XARS, which differ from the PHFIT²³ cross sections used in SKIRT, REFLE X, and XSPEC. The XARS cross sections are 10 to 15% higher for photon energies below 12.5 keV, causing the corresponding soft X-ray fluxes to be lower for obscured sightlines. The induced flux reduction is expected to further depend on the column density and the photon energy, as shown in the bottom right panel of Fig. 17.

At higher photon energies, the total gas extinction is fully dominated by Compton scattering (see Figs. 4 and 7), which is treated consistently in both codes and therefore reconciles the SKIRT and XARS results. In particular, the Compton reflection hump peaking at about 40 keV is nicely recovered on top of the unobscured input spectrum with both codes, as shown in the left top panel of Fig. 17. This same figure shows a clear S K_α line in the SKIRT results at 2.3 keV, which is the strongest fluorescent line that is not implemented in XARS, arguably more important than the available Ne and Mg transitions. Finally, we find an excellent agreement between both codes in the 6.0 to 7.5 keV spectral region for heavily obscured sightlines, reproducing the characteristic spectral shape of heavily obscured AGN (top right panel of Fig. 17, for N_H = 10²⁴ cm⁻²). Overall, less Monte Carlo noise is observed in the SKIRT results, which can be linked to the implemented optimisation techniques (Sect. 2).

Fig. 17 Seventh SKIRT benchmark (Sect. 4.2.4), comparing SKIRT simulations of a torus model against the corresponding XARS spectra. The opening angle of the wedge torus is 60°. The equatorial hydrogen column density NH of the torus is indicated on the figure, for each simulation. Left: total photon flux, for a line of sight not intersecting the obscuring torus (i = 45°). Right: total photon flux, looking through the obscuring torus (i = 75°). Both codes are in excellent agreement. The small offset below 12.5 keV is explained by differences in the atomic data (see text).

5 Demonstration

We illustrate the 3D capabilities of SKIRT in the X-ray regime by simulating a set of clumpy torus models, where a variable fraction of the torus material is concentrated in clumps. These clumps are randomly distributed over the torus, forming a truly 3D transfer medium with no spatial symmetries to reduce the dimensionality of the radiative transfer simulations.

The baseline smooth torus model is a uniform-density wedge torus of cold gas with an opening angle of 60^°, and an equatorial hydrogen column density of 5 × 10²² cm⁻², which is centred around an X-ray point source with г = 1.8, E_cut = 200 keV, and L_2–10 = 2.8 × 10⁴² erg s⁻¹ (see Sect. 4). The radial extent of the torus is 15 pc, and the gas abundances are assumed to be solar, with bound-electron scattering enabled. We then consider the corresponding two-phase torus model by allocating 50% of the torus mass to 1000 spherical clumps of 0.5 pc radius, and the clumpy torus model that has all its mass in clumps (see also Stalevski et al. 2012). Within each clump, we assume a cubic spline radial density profile. For details on the implementation of such a clumpy geometry we refer to Baes & Camps (2015).

Figure 18 shows the SKIRT X-ray images for the smooth, two-phase, and clumpy torus models described above. The smooth torus model (Fig. 18, left) exhibits a colour asymmetry between the illuminated back and front side of the torus, which directly relates to the X-ray scattering physics in the gas (Sect. 3.2.3). The two-phase and clumpy torus models show more complex morphologies, and demonstrate how individual clumps can be properly resolved in 3D radiative transfer simulations with SKIRT. These clumps form dense scattering sites inside the torus, which are not easily penetrated by soft X-ray photons, producing radially extended shadow lanes (Fig. 18, middle). For the two-phase and clumpy torus models, the column densities of individual clumps reach 9.0 × 10²² cm⁻² and 1.8 × 10²³ cm⁻², respectively. Overall, these images illustrate how the entire torus structure contributes to the X-ray reflection spectrum, through multiple scattering inside the 3D transfer medium.

The corresponding X-ray spectra are shown in Fig. 19, demonstrating how clumpy torus models allow some soft X-ray photons to escape through the less dense inter-clump material. We note how the simulated differences between these spectra might be difficult to distinguish in real observations, due to the presence of additional strong components in the soft X-ray range, such as Thomson scattered emission originating from larger distances (see Gupta et al. 2021, and Fig. 19). Finally, we illustrate how torus models lose their axial symmetry when making them clumpy, by observing one clumpy torus (f_clumps = 100%) from different azimuthal viewing angles in Fig. 20. Depending on the distribution of clumps along the line of sight to the observer, the observed X-ray spectra can differ drastically. For example, the line of sight corresponding to ∅ = 80° does not intersect any clumps, whereas the one corresponding to ∅ = 10° intersects several clumps, which leads to an attenuation of the soft X-ray flux at 1 keV by more than 3 orders of magnitude. We refer to future SKIRT work for a more in-depth systematic study of different 3D torus geometries, under various physical assumptions.

Fig. 18 SKIRT X-ray images of a smooth torus model (left), two-phase torus model (middle), and clumpy torus model (right), observed at i = 75°. The RGB colours correspond to the integrated fluxes in the given ranges. Clumpy AGN tori represent truly 3D transfer geometries with no spatial symmetries.

Fig. 19 SKIRT X-ray spectra for torus models with different clump mass fractions fclumps The red, blue, and black spectra correspond to the smooth, two-phase, and clumpy tori in Fig. 18, respectively. The grey dotted line represents additional Thomson scattered emission at 1% of the intrinsic continuum (see Gupta et al. 2021).

Fig. 20 SKIRT X-ray spectra for a clumpy torus model (fclumps = 100%), observed from different azimuthal viewing angles (with i = 75◦ fixed).

6 Discussion and conclusions

We presented the new X-ray functionalities of the MCRT code SKIRT, for studying cold dusty gas in complex 3D transfer media. Focussing on continuum radiative transfer in the 0.1 to 500 keV range, we introduced Compton scattering on free electrons (Sect. 3.1), photo-absorption (Sect. 3.2.1) and fluorescence (Sect. 3.2.2) by cold atomic gas, scattering on bound electrons (Sect. 3.2.3), and extinction by dust (Sect. 3.3.1). This includes a novel treatment of extreme-forward scattering by dust (Sect. 3.3.2), and a detailed description of anomalous Rayleigh scattering (Sect. 3.2.3). By extending the SKIRT code with these X-ray processes, we obtained an X-ray radiative transfer code with all features of the established SKIRT framework (Sect. 2), fully optimised to operate in general 3D geometries. Furthermore, by extending the wavelength domain of SKIRT into the X-ray range, we enabled radiative transfer simulations that can self-consistently cover the X-ray band plus the IR-to-UV wavelength range, linking X-ray reprocessing to dust modelling.

The main motivation for extending the advanced treatment of absorption, scattering, and re-emission into the X-ray range, was to study the structure of AGN circumnuclear media based on their reprocessed X-ray emission. With SKIRT, we can now experiment with X-ray obscuration models of arbitrary complexity, to challenge the classical torus paradigm in a computation-efficient way. Indeed, the geometrical flexibility of SKIRT allows for testing a broad suite of spectral models featuring clumps and filaments, polar structures, and distributions imported from hydrodynamical simulations. Next to this 3D-efficiency, SKIRT offers a complete set of X-ray physics for modelling gas and dust in AGN circumnuclear media, exceeding what is currently implemented in most X-ray spectral models. Apart from AGN, we foresee additional applications for the SKIRT code in the X-ray band, for example in modelling galactic sources such as novae and X-ray binaries, or modelling the X-ray SEDs of galaxies. Finally, the advanced treatment of extreme-forward scattering by dust makes SKIRT a powerful tool for studying X-ray scattering halos formed in complex dust geometries.

We verified our X-ray implementation by reproducing some well-established XSPEC models, and found that our SKIRT results were perfectly consistent with the standard spectral models for free-electron extinction (CABS, Sect. 4.1.1), cold-gas absorption (PHABS, Sect. 4.1.1), fluorescence (PEXMON, Sect. 4.1.2), Compton scattering (PEXMON, Sect. 4.1.2), and dust extinction (ISMDUST, Sect. 4.1.3). These 1D XSPEC models incorporate few non-trivial radiative transfer effects, but assure a correct implementation of the basic X-ray processes in free-electron, cold-gas, and dust media.

We performed the first dedicated benchmark of X-ray torus models, comparing five popular X-ray radiative transfer codes (MYTORUS, REFLEX, BORUS, XARS, and SKIRT). We found an excellent agreement between the SKIRT results and the MYTORUS (Sect. 4.2.1) and REFLEX (Sect. 4.2.2) model spectra, for various parameter combinations. REFLEX implements bound-electron scattering just as SKIRT and uses the same atomic data, which makes it ideal for benchmarking the new X-ray processes in SKIRT. We recovered both the reflected X-ray continuum (Fig. 14), the fluorescent line fluxes (Table 2), and the Fe K_α Compton shoulder (Fig. 15), forming a convincing validation of the underlying X-ray implementation.

We found a good agreement between the SKIRT and XARS results modelling a wedge torus medium, which is one of the many transfer geometries that was released with XARS (Sect. 4.2.4). Compared to SKIRT, we observed slightly more photo-absorption in the XARS spectra at soft X-ray energies, which we traced back to a difference in the adopted photo-absorption cross sections. At higher energies, we found an excellent agreement between the Compton reflection humps as obtained with both codes (Fig. 17), reconciling both implementations of free-electron scattering. In Sect. 4.2.3, we found significant discrepancies between the SKIRT and BORUS simulation results, which were partly understood. These differences, observed at low and high energies (Fig. 16), plus the minor soft X-ray offset relative to XARS, illustrate the complexity of X-ray radiative transfer and motivate the need for a robust framework that can handle non-linear radiative transfer effects in full 3D.

Most models featured in our X-ray torus benchmark showed more numerical noise than the corresponding SKIRT results, which should be linked to the MCRT optimisations implemented in SKIRT (Sect. 2). In Monte Carlo simulations, the signal-to-noise ratio is expected to scale as the square root of the total number of photon packets, while increasing the total number of photons affects the simulation time linearly. MCRT acceleration techniques allow for reaching similar noise levels with a reduced number of photons, thus accelerating the simulation process. For a fixed number of photons, this means that SKIRT can run radiative transfer simulations at a higher spectral resolution, while maintaining reasonable noise levels over more spectral bins. in this work, all SKIRT spectra were calculated at a higher spectral resolution than the corresponding reference model, and contain significantly less noise.

In addition to this MCRT-efficiency, the SKIRT code is highly optimised in a purely computational sense, resulting in shorter simulation times for a fixed number of photon packets. This reduces the computational cost of radiative transfer simulations, which is important for calculating large model libraries. in this work, each SKIRT simulation was run on a 2.3 GHz 8-core laptop in just a few hours, which is remarkably short considering that each simulation featured as many as 5 × 10⁹ photons. Furthermore, these 5 × 10⁹ photons are directly contributing to the simulation output, which should not be compared to the total number of photons escaping the model in all directions. increasing the total number of photons linearly affects the simulation time, but does not affect the memory requirements (assuming the same host computing system). increasing the spectral resolution has a very limited effect on the overall memory usage, but does require launching extra photons to properly sample each energy bin. Finally, SKIRT scales efficiently on multiple cores and multiple nodes, allowing to significantly increase the computational resources when needed (Verstocken et al. 2017).

MCRT simulations are mostly limited by extreme optical depths of τ > 20, depending on the model geometry and the relative placement of sources and detectors (Camps & Baes 2018). For hard X-ray photons, this corresponds to column densities of N_H > 2 × 10²⁵ cm⁻² (see Fig. 7), which is sufficiently high to cover the entire column density distribution of AGN in the local universe (Ricci et al. 2015). Indeed, SKIRT can produce transmission spectra for line-of-sight column densities of N_H = 10²⁵ cm⁻² without significant numerical noise, as shown in Fig. 10. We note how these density constraints hold for the transmitted flux component only, with no similar restrictions for the reflected flux, as demonstrated in Sect. 4.1.2 where reflection on an dense slab of N_H = 10²⁷ cm⁻² was modelled.

In Sect. 5, we demonstrated the full-3D capabilities of SKIRT in the X-ray regime by simulating X-ray images and spectra of clumpy torus models, illustrating how SKIRT will be used to build spectral models based on complex transfer geometries. All configuration elements (such as medium geometries) needed to run these clumpy torus models, plus the benchmark simulations presented in Sect. 4, are built into the public master branch of the SKIRT code, and can used without any further modifications.

The aforementioned features make SKIRT a promising new tool for modelling the complex structure of AGN circumnuclear media based on their X-ray emission. With SKIRT, we can predict fluorescent lines and Compton shoulder shapes that cannot be observed with the current generation of X-ray observatories, but should appear in upcoming microcalorimeter observations (Barret et al. 2018; Tashiro et al. 2020). Arguably, the access to full-3D models of complex obscuring structures will be crucial for modelling high-resolution XRISM/Resolve data, where we expect informative spectral features that cannot be reproduced with 1D models. The presented update of the SKIRT code allows for uncomplicated access to a wide suite of 3D X-ray models for obscured AGN that can easily be tested and modified, with limited computational efforts.

Acknowledgements

We wish to thank the anonymous referee for their careful reading and constructive comments that helped improve the presentation of this work. B.V. acknowledges support by the Fund for Scientific Research Flanders (FWO-Vlaanderen, project 11H2121N). M.S. acknowledges support by the Science Fund of the Republic of Serbia, PROMIS 6060916, BOWIE and by the Ministry of Education, Science and Technological Development of the Republic of Serbia through the contract no. 451-03-9/2022-14/200002. We wish to thank K. A. Arnaud, B. T. Draine, A. A. Zdziarski, L. Corrales and M. Balokovic for helpful discussions.

References

All Tables

All Figures

	Fig. 1 Total photo-absorption cross section for a cold gas medium with Anders & Grevesse (1989) abundances and no dust. The K- and L-shells dominate the absorption over the entire X-ray range. The prominent absorption edge structure reflects the atomic composition of the gas.
In the text

	Fig. 2 Main contributions to the total photo-absorption cross section of a cold gas (see Fig. 1), assuming Anders & Grevesse (1989) abundances and no dust.
In the text

	Fig. 3 Fluorescent yields for all fluorescent transitions that are implemented in SKIRT, as provided by Perkins et al. (1991). See Table 1 for the transition definitions.
In the text

	Fig. 4 Bound-electron scattering cross section of a cold gas medium, assuming Anders & Grevesse (1989) abundances and no dust. The 1.21 free-electron approximation is shown in grey for comparison (see text).
In the text

	Fig. 5 Normalised scattering phase functions for bound-electron scattering in cold gas, at a photon energy of 3.3 keV. At this photon energy, bound-Compton scattering and Rayleigh scattering are equally important, see Fig. 4. The free-electron phase function at 3.3 keV is shown in grey for comparison.
In the text

	Fig. 6 Anomalous Rayleigh scattering cross sections of H, Be, Mg, and Zn. The corresponding smooth Rayleigh cross sections (see Eq. (18)) are shown in black.
In the text

	Fig. 7 Total gas extinction cross section for Anders & Grevesse (1989) abundances and no dust. The effective photo-absorption (blue), bound-electron scattering (green), and effective fluorescent scattering (yellow) contributions are explicitly shown, to illustrate how these processes shape the X-ray spectrum in radiative transfer simulations.
In the text

	Fig. 8 Dust extinction cross sections in the X-ray range, as implemented in SKIRT. The atomic absorption edges of the elements that make up the dust grains are indicated in grey.
In the text

	Fig. 9 Dampened relative bias factors 〈b〉 for dealing with extreme forward dust scattering. For δ = 4°, we maintain 〈b〉 ≲ 820.
In the text

	Fig. 10 First SKIRT benchmark (Sect. 4.1.1), comparing three SKIRT simulations with the PHABS and CABS models for cold-gas and free-electron extinction. The hydrogen column density NH is indicated on the figure, while the electron column density isNe = 1.21 × NH (see text).
In the text

	Fig. 11 Second SKIRT benchmark (Sect. 4.1.2), comparing three SKIRT simulations with the PEXMON model for reflection by neutral material. The inclination relative to the reflecting plane is indicated on the figure. The 10% discrepancy in the 12–60 keV range was expected (see text).
In the text

	Fig. 12 Third SKIRT benchmark (Sect. 4.1.3), comparing SKIRT with the ISMDUST model for dust extinction. The dust mass column density Σd is indicated on the figure. Both models are in good agreement above the C K-edge at 0.280 keV (see text).
In the text

Fig. 13 Fourth SKIRT benchmark (Sect. 4.2.1), comparing SKIRT simulations of a torus model against the corresponding MYTORUS model spectra. The opening angle of the ring torus is 60°. The equatorial hydrogen column density NH of the torus is indicated on the figure, for each simulation. Left: scattered photon flux, for a line of sight not intersecting the obscuring torus (i = 45°). Right: scattered photon flux, looking through the obscuring torus (i = 75°). We find an excellent agreement between both codes, for all column densities and sightlines.

In the text

Fig. 14 Fifth SKIRT benchmark (Sect. 4.2.2), comparing SKIRT simulations of a torus model against the corresponding RXTORUS model spectra as calculated with REFLEX. The opening angle of the ring torus is 60°. The equatorial hydrogen column density NH of the torus is indicated on the figure, for each simulation. Left: scattered photon flux, for a line of sight not intersecting the obscuring torus (i = 45°). Right: scattered photon flux, looking through the obscuring torus (i = 75°). We find an excellent agreement between both codes, for all column densities and sightlines.

In the text

	Fig. 15 Zoom on the 6.0 to 7.5 keV reflection spectrum of a ring torus model with NH = 1024 cm−2 (same model as Fig. 14, left). Table 2 lists the integrated line fluxes for the three most prominent fluorescent lines. We note how SKIRT can spectrally resolve the Kα1 and Kα2 sublines.
In the text

Fig. 16 Sixth SKIRT benchmark (Sect. 4.2.3), comparing SKIRT simulations of a torus model against the corresponding BORUS02 model spectra as calculated with BORUS. The opening angle of the wedge torus is 60°. The equatorial hydrogen column density NH of the torus is indicated on the figure, for each simulation. Left: scattered photon flux, for a line of sight not intersecting the obscuring torus (cos i = 0.75). Right: scattered photon flux, looking through the obscuring torus (cos i = 0.25). Above 30 keV, we find a significant mismatch between both codes for all sightlines. For unobscured sightlines (cos i = 0.75), the simulation results also diverge for energies below 10 keV.

In the text

Fig. 17 Seventh SKIRT benchmark (Sect. 4.2.4), comparing SKIRT simulations of a torus model against the corresponding XARS spectra. The opening angle of the wedge torus is 60°. The equatorial hydrogen column density NH of the torus is indicated on the figure, for each simulation. Left: total photon flux, for a line of sight not intersecting the obscuring torus (i = 45°). Right: total photon flux, looking through the obscuring torus (i = 75°). Both codes are in excellent agreement. The small offset below 12.5 keV is explained by differences in the atomic data (see text).

In the text

	Fig. 18 SKIRT X-ray images of a smooth torus model (left), two-phase torus model (middle), and clumpy torus model (right), observed at i = 75°. The RGB colours correspond to the integrated fluxes in the given ranges. Clumpy AGN tori represent truly 3D transfer geometries with no spatial symmetries.
In the text

	Fig. 19 SKIRT X-ray spectra for torus models with different clump mass fractions fclumps The red, blue, and black spectra correspond to the smooth, two-phase, and clumpy tori in Fig. 18, respectively. The grey dotted line represents additional Thomson scattered emission at 1% of the intrinsic continuum (see Gupta et al. 2021).
In the text

	Fig. 20 SKIRT X-ray spectra for a clumpy torus model (fclumps = 100%), observed from different azimuthal viewing angles (with i = 75◦ fixed).
In the text