X-ray polarisation in AGN circumnuclear media. Polarisation framework and 2D torus models
Bert Vander Meulen Peter Camps Djordje Savic Maarten Baes Giorgio Matt Marko Stalevski
 


   

   Polarisation framework and 2D torus models

   X-ray polarisation in AGN circumnuclear media

   B. Vander Meulen et al.

   
   Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281 S9, 9000 Gent, Belgium

              bert.vandermeulen@ugent.be
              Institut d’Astrophysique et de Géophysique, Université de Liège, Allée du 6 Août 19c, 4000 Liège, Belgium
              Astronomical Observatory, Volgina 7, 11060 Belgrade, Serbia
              Dipartimento di Matematica e Fisica, Universitá degli Studi Roma Tre, via della Vasca Navale 84, 00146 Roma, Italy
        
   




  
  
   Cold gas and dust reprocess the central X-ray emission of active galactic nuclei (AGN), producing characteristic spectro-polarimetric features in the X-ray band. The recent launch of IXPE allows for observations of this X-ray polarisation signal, which encodes unique information on the parsec-scale circumnuclear medium of obscured AGN. However, the models for interpreting these polarimetric data are under-explored and do not reach the same level of sophistication as the corresponding spectral models.
  
   We aim at closing the gap between the spectral and spectro-polarimetric modelling of AGN circumnuclear media in the X-ray band by providing the tools for simulating X-ray polarisation in complex geometries of cold gas alongside X-ray spectra.
  
   We lay out the framework for X-ray polarisation in 3D radiative transfer simulations and provide an implementation to the 3D radiative transfer code SKIRT, focussing on (de)polarisation due to scattering and fluorescent re-emission. As an application, we explored the spectro-polarimetric properties of a 2D toroidal reprocessor of cold gas, modelling the circumnuclear medium of AGN.
  
   For the 2D torus model, we find a complex behaviour of the polarisation angle with photon energy, which we interpret as a balance between the reprocessed photon flux originating from different sky regions, with a direct link to the torus geometry. We calculated a large grid of AGN torus models and demonstrated how spatially resolved X-ray polarisation maps could form a useful tool for interpreting the geometrical information that is encoded in IXPE observations. With this work, we release high-resolution AGN torus templates that simultaneously describe X-ray spectra and spectro-polarimetry for observational data fitting with XSPEC.
  
   The SKIRT code can now model X-ray polarisation simultaneously with X-ray spectra and provide synthetic spectro-polarimetric observations for complex 3D circumnuclear media, with all features of the established SKIRT framework available.

    

   X-ray polarisation in AGN circumnuclear media
    Bert Vander Meulen
          1
          Peter Camps
          1
          Đorđe Savić
          2, 3
          Maarten Baes
          1
          Giorgio Matt
          4
          Marko Stalevski
          1, 3
          
    Received May 17, 2024; accepted July 29, 2024
===================================================================================================================================================================================================================================







§ INTRODUCTION



Active galactic nuclei (AGN) are the compact central regions of massive galaxies whose excessive brightness across the electromagnetic spectrum is powered by the accretion of gas and dust onto a supermassive black hole (SMBH) <cit.>. AGN are the most luminous persistent sources in the Universe and play a crucial role in galaxy evolution by pumping energy and momentum into the interstellar gas and launching powerful jets and outflows. AGN feedback is one of the most important and least understood ingredients of galaxy evolution theories, with several observed correlations indicating that AGN and their host galaxy co-evolve and regulate each other’s growth <cit.>.


According to the unified AGN structure model, the observed dichotomy in AGN types is explained by a large-scale toroidal structure of gas and dust in the equatorial plane, which causes line-of-sight obscuration, depending on the observer's viewing angle <cit.>. This `torus' of gas and dust is then responsible for the extinction at optical and UV wavelengths, the thermal dust re-emission in the infrared, and reprocessing in the X-ray band. Furthermore, it explains the absence of broad optical lines in obscured AGN and the appearance of these lines in polarised light <cit.>. This circumnuclear torus medium could further be important as an accretion reservoir fuelling the active SMBH or as a direct probe on AGN feedback <cit.>.


Lately, a new picture has emerged of the dust structure in local active galaxies, challenging the classical `dusty torus' paradigm. Spectral modelling suggests that the circumnuclear medium has a more complex three-dimensional structure with clumps and filaments <cit.>, which is further supported by observations of AGN variability <cit.>. Furthermore, high-angular resolution imaging in the mid-infrared has shown that the circumnuclear medium is also extended in the polar direction, as opposed to a purely equatorial torus  <cit.>.


The AGN circumnuclear medium could be further explored through X-ray observations, as most AGN spectra show a strong X-ray component, which is produced by Compton up-scattering in a corona of hot electrons close to the SMBH <cit.>. These coronal X-rays are then reprocessed by the circumnuclear material, producing characteristic spectral features that form a powerful probe regarding the distribution of gas and dust in local AGN <cit.>. X-rays have a high penetrating power, and therefore they shed light on the most obscured episodes of SMBH accretion. Indeed, a large population of obscured AGN is only revealed through X-ray observations, with spectra that are shaped by reprocessing in the AGN torus <cit.>. As many of these sources cannot be spatially resolved, they contribute to the cosmic X-ray background <cit.>.


The most prominent features of X-ray reprocessing by the AGN torus are the narrow Fe Kα line at 6.4 keV <cit.> and the Compton reflection hump peaking at about 30 keV <cit.>. These features directly probe the cold gas and dust surrounding AGN and are commonly used to constrain the geometry of the reprocessing medium <cit.>. Indeed, recent radiative transfer studies have demonstrated that X-ray spectra of obscured AGN carry detailed information on the 3D distribution of obscuring gas and dust, tracing clumpy structures <cit.> and polar extended material <cit.> in the circumnuclear medium, even outside of the line of sight.


Two additional observables can be extracted from the X-ray emission emerging from obscured AGN when dedicated polarisation instrumentation is available. The polarisation angle and polarisation degree encode complementary information on the circumnuclear medium, which can be used to constrain the geometry of the torus and its orientation relative to the host galaxy. The recent launch of the Imaging X-ray Polarimetry Explorer <cit.> introduced X-ray polarimetry as a new tool to study AGN in the 2-8 keV band, with five radio-quiet AGN that have been observed in the last two years: MCG-05-23-16 <cit.>, the Circinus galaxy <cit.>, NGC4151 <cit.>, IC4329A <cit.>, and NGC1068 <cit.>. The next generation of X-ray polarisation missions, with the X-ray Polarimeter Satellite (XPoSat) <cit.> and the enhanced X-ray Timing and Polarization mission (eXTP) <cit.>, further indicate a promising future for X-ray polarimetry.


The recent developments in observational X-ray polarimetry motivate the need for more advanced polarisation models that are based on 3D radiative transfer simulations. However, spectro-polarimetric X-ray models for the AGN torus are under-explored and do not reach the same level of sophistication as the corresponding X-ray spectral models <cit.>.


Historically, one of the first X-ray polarisation models was the wedge torus model by <cit.>, which has been used for the interpretation of the 772 ks of IXPE data on Circinus AGN <cit.> and for the exploration of the binary geometry of GRS 1915+105 <cit.>. However, this model allows for little geometrical flexibility.


The STOKES code <cit.> on the other hand offers a more flexible radiative transfer framework, and it has been used to make polarisation predictions for a range of 3D torus geometries <cit.>. STOKES applies a Monte Carlo technique to model scattering-induced polarisation in the X-ray, UV, optical, and infrared bands. Recently, the STOKES code has been used to model X-ray polarisation in a parsec-scale equatorial torus, assuming neutral and partially ionised gas <cit.>, which has been applied to the 1.15 Ms IXPE observation of NGC1068 <cit.>.


Most recently, the MONACO code <cit.> has been used to predict the X-ray polarisation signal of Circinus AGN, which was compared to the IXPE observational data <cit.>. MONACO, which builds on the Geant4 simulation toolkit <cit.>, implements X-ray polarisation in cold neutral material, which has been applied to post-process the 3D hydrodynamical torus simulations by <cit.> modelling the circumnuclear gas in the Circinus galaxy.


With this work, we aim to close the gap between the X-ray spectral modelling and X-ray spectro-polarimetric modelling of AGN circumnuclear media by setting up 3D radiative transfer simulations that simultaneously predict X-ray polarisation and X-ray spectra. For this, we focus on the well-established SKIRT code <cit.> so that our simulations are highly efficient in terms of computational runtime, allowing complex 3D models to be explored in a short time. Furthermore, the SKIRT code offers an unmatched geometrical flexibility for setting up radiative transfer simulations in full 3D, which has now become available to X-ray polarisation modelling. This work introduces the necessary tools for modelling X-ray polarisation in a general 3D context and presents an implementation to the SKIRT code, which is publicly available. As a first application, we explored the spectro-polarimetric properties of a simple 2D torus model. Future work will focus on more complex 3D geometries beyond the classical AGN torus.


The goal of this work is to lay out the X-ray polarisation framework for 3D radiative transfer codes and investigate the X-ray polarisation observables of a classical AGN torus model. The outline for this paper is as follows: In Sect. <ref>, we introduce the polarisation framework. In Sect. <ref>, we provide an implementation to the 3D radiative transfer code SKIRT. In Sect. <ref>, we study the polarisation observables of a 2D torus model and present the corresponding torus templates for observational data fitting. We discuss our findings in Sect. <ref> and summarise our results in Sect. <ref>.



§ X-RAY POLARISATION FRAMEWORK




 §.§ Stokes vector formalism


In Monte Carlo radiative transfer simulations, the polarisation state of the (discretised) radiation field is most conveniently characterised in terms of the four Stokes parameters <cit.>. Together, these parameters form a Stokes vector , which is traced and updated throughout the simulated transfer medium for each `photon packet':

    = [ I; Q; U; V ].

Stokes parameter I represents the total intensity of the photon packet, while Q and U describe linear polarisation, and V describes circular polarisation. The four Stokes parameters encode all polarisation information of the radiation field (except for its phase), so that the linear polarisation degree P_L can be reconstructed as

    P_L = √(Q^2 + U^2)/I,

and the linear polarisation angle γ can be found as

    γ = 1/2arctan_2(U/Q),

with arctan_2 being the inverse tangent function that preserves the quadrant, also noted as atan2(U, Q).

Stokes vectors  are defined relative to a reference direction  that is perpendicular to the photon propagation direction . Positive Q-values describe linear polarisation in the reference direction, while negative Q-values describe linear polarisation perpendicular to this reference direction. Equivalently, Stokes U describes linear polarisation along orthogonal axes rotated over  45^∘ with respect to the reference direction. One has


    Q    = I P_Lcos2γ,
    
    U    = I P_Lsin2γ.


For the remainder of this work, we focus on linear polarisation, as IXPE observations do not capture the Stokes V-component. Nevertheless, our polarisation framework includes a functional implementation of circular polarisation (which is not discussed).

When the reference direction  is rotated by an angle φ about  to a new reference direction _rot, the Stokes parameters transform as

    _rot = R(φ) ,

with R(φ) being the rotation matrix, which is given as

    R(φ) =
    [      1      0      0      0;      0  cos2φ  sin2φ      0;      0 -sin2φ  cos2φ      0;      0      0      0      1 ],

mixing the Q and U parameters. Combining Eq. (<ref>), (<ref>) and (<ref>), we obtain


    Q_rot   = I P_Lcos2(γ - φ),
    
    U_rot   = I P_Lsin2(γ - φ),


so that we have P_L_rot=P_L and γ_rot=γ - φ. This rotation transformation is important to describe scattering interactions relative to a reference direction that lays in the scattering plane and to record Stokes vectors relative to the north direction of the observer frame (following the IAU conventions).


In polarised radiative transfer, photon interactions generally depend on the polarisation state, meaning that the processes are described in terms of the Stokes parameters of that photon. In addition, scattering can modify the polarisation state of a photon, which can be described as a matrix multiplication on  with the corresponding Müller scattering matrix M. In this work, Müller matrices M(θ, x) are a function of the scattering angle θ and the incoming photon energy x = E/m_ec^2. Considering that the reference direction of the incoming photon should first be rotated to an intermediate reference direction _rot in the scattering plane, the Stokes vector ' after scattering can be obtained as

    '=M(θ,x) R(φ) ,

where φ is the angle by which the initial reference direction  must be rotated about  to end up in the scattering plane as _rot (see Fig. <ref>). Finally, we note that Stokes vector ' refers to a new reference direction ', which is just the intermediate reference direction _rot, rotated over the scattering angle θ in the scattering plane (i.e. a rotation about the scattering plane normal ×_rot), to assure that the new reference direction ' is perpendicular to the new propagation direction '. The 3D geometry of a scattering interaction is shown in Fig. <ref>, with unit direction vectors, unit reference vectors, and rotation angles indicated. For more details, we refer to <cit.>. Applying Rodrigues' rotation formula, we obtain an expression for ' as

    '=cosφcosθ  + sinφcosθ(×) -sinθ .


Various normalisation conventions are in use for the Müller matrices M (e.g. different conventions can be found in <cit.>, <cit.> and <cit.>). However, polarisation observables such as P_L and γ do not depend on this absolute normalisation factor, as they are calculated from Stokes parameter ratios, causing the normalisation factor to cancel out. The absolute polarised fluxes Q and U can be recovered from the normalised total flux I (similar to Eq. (<ref>)), which is discussed in Sect. <ref>.



 §.§ X-ray radiative processes


In this work, we focus on radiative transfer simulations in cold atomic gas, modelling the circumnuclear medium that is causing most of the X-ray extinction in obscured AGN. Furthermore, this material is responsible for the distinct X-ray reflection features observed in Compton-thick AGN, such as the narrow Fe Kα line at 6.4 keV and the Compton reflection hump at about 30 keV. In particular, we consider photo-absorption by neutral atomic gas with self-consistent fluorescent line emission and bound-electron scattering.

Photo-absorption does not depend on the polarisation state of the incoming photon. Therefore, it can be implemented using the standard <cit.> cross-sections and a custom abundance table. Subsequently, the absorbed photon energy can be re-emitted as a fluorescent line photon, with a probability given by the fluorescent yield of that atom <cit.>. These line photons are unpolarised, regardless of the initial polarisation state of the absorbed photon <cit.>. Fluorescence thus resets the Stokes vector, and the sequence of photo-absorption followed by fluorescent re-emission effectively acts as a depolarisation process.

In cold atomic gas, X-ray photons are scattered by the electrons that are bound to the neutral gas atoms, which can be described in terms of atomic form factors and incoherent scattering functions <cit.>. However, this bound-electron scattering can be reasonably well approximated as Compton scattering on Z free electrons per neutral gas atom, which we assume in this work. For <cit.> solar gas abundances, this corresponds to 1.21 free electrons per H-atom, which changes for other abundance tables. This free-electron approximation, applied to X-ray polarisation in a general 3D context, is the focus of this work, and we refer to future work for a treatment of polarised bound-electron scattering <cit.>.

In the next section, we describe the details of polarised Compton scattering (approximating bound-electron scattering), which is polarisation dependent and updates the polarisation state of the incoming photon. In particular, we focus on the total scattering cross-section and the scattering phase function for polarised photons and the Müller matrix for updating the Stokes parameters after a scattering interaction. This framework could also be used to model media of free electrons only, in addition to the cold gas that is the focus here.



 §.§ Polarised Compton scattering


Compton scattering describes the inelastic scattering of high-energy photons on free electrons, which is often used to approximate bound-electron scattering in cold atomic gas <cit.>. The Müller matrix for Compton scattering is given by <cit.>, using the same sign convention as the IAU <cit.>:

    M(θ, x) ∝[ S_11(θ, x) S_12(θ, x)          0          0; S_12(θ, x) S_22(θ, x)          0          0;          0          0 S_33(θ, x)          0;          0          0          0 S_44(θ, x) ],

with θ being the scattering angle and x = E/m_ec^2 the photon energy scaled to the electron rest energy. Moreover, we assume an isotropic distribution for the electron spin directions, which causes the non-diagonal matrix elements in the fourth row and fourth column to be zero <cit.>. The non-zero matrix elements of M(θ, x) are



    S_11(θ, x)    = C^3(θ, x) + C(θ, x) -C^2(θ, x)sin^2θ,
    S_12(θ, x)    = -C^2(θ, x)sin^2θ,
    S_22(θ, x)    = C^2(θ, x) (1+cos^2θ),
    S_33(θ, x)    = 2C^2(θ, x)cosθ,
    S_44(θ, x)    = (C^3(θ, x) + C(θ, x)) cosθ,


with C(θ, x) being the Compton factor:

    C(θ, x) = (1+x  (1-cosθ))^-1.

In the low energy limit, C(θ, x)→ 1, so that Eq. (<ref>) converges to the Müller matrix for non-relativistic Thomson scattering.

Matrix element S_11(θ, x) relates the scattered intensity I' to the incoming intensity I in case of unpolarised photons and is proportional to the <cit.> differential cross-section for unpolarised Compton scattering. For (partially) polarised radiation, S_12(θ, x) introduces an azimuthal modulation to the Klein-Nishina formula, which depends explicitly on the polarisation state of the incoming photon. Together, S_11(θ, x) and S_12(θ, x) make up the differential cross-section for polarised Compton scattering:

    dσ/dΩ(θ, φ, x, ) = 3σ_T/16π[S_11(θ, x) + S_12(θ, x) P_Lcos2(φ-γ)],

with σ_T≈ 6.65 × 10^-25 cm^2 being the Thomson cross-section and φ the azimuthal scattering angle relative to the polarisation reference direction  (so that φ-γ is the azimuthal scattering angle relative to the incoming polarisation direction ).

For fully polarised photons (P_L=1), Eq. (<ref>) reduces to

    dσ/dΩ(θ, φ, x, )|_P_L=1 = 3 σ_T/16π[ C^3 + C -2 C^2sin^2θcos^2(φ-γ)],

with C ≡ C(θ, x). In the IXPE range, the angular dependence of C(θ, x) is weak, and Eq. (<ref>) is roughly symmetric around the incoming polarisation direction , as the factor sinθ cos(φ-γ) is just the cosine between ' and . For a full derivation from quantum electrodynamics, we refer to <cit.>.

By integrating Eq. (<ref>) over the unit sphere, we obtain the total cross-section for polarised Compton scattering as

    σ(x) = 3 σ_T/4[1+x/(1+2x)^2 + 2/x^2 + (1/2x-x+1/x^3) ln(2x+1) ],

which is the exact same expression as the total cross-section for unpolarised Compton scattering. This cross-section does not depend on the polarisation state, meaning that polarisation effects do not influence the Compton scattering efficiency. Furthermore, this cross-section is roughly constant over the IXPE energy range, decreasing from 0.99 σ_T at 2 keV to 0.97 σ_T at 8 keV.


Normalising Eq. (<ref>) on the unit sphere, we obtain the phase function for polarised Compton scattering, shown in Fig. <ref>:

    Φ(θ, φ, x, )    = 3σ_T/16πσ(x)[S_11(θ, x) + S_12(θ, x) P_Lcos2(φ-γ)]
       = Ψ(θ, x) - 3σ_T  C^2(θ, x)/16πσ(x)sin^2θ P_Lcos2(φ-γ).

This phase function Φ(θ, φ, x, ) is the sum of the standard phase function for unpolarised Compton scattering Ψ(θ, x)∝S_11(θ, x) and an azimuthal modulation that is characteristic for polarised Compton scattering. The azimuthal cosine modulation introduces a bias towards scattering into the plane perpendicular to the polarisation direction (φ= γ + 90^∘ and φ= γ + 270^∘), with an amplitude that is roughly proportional to sin^2θ× P_L. Therefore, the modulation is most pronounced at θ=90^∘, where it maximally reduces the probability for scattering into the polarisation direction (φ= γ or φ= γ + 180^∘). As expected, the modulation strength is linearly proportional to the polarisation degree P_L, as only polarised photons experience the modulation effect.

The left panel of Fig. <ref> shows the Compton phase function at 2 keV (solid line) and 8 keV (dashed line), illustrating how the azimuthal modulation does not depend on the photon energy in the IXPE range. The only difference is the slightly stronger preference for forward scattering at higher energies, which is regulated by the polarisation-independent phase function term Ψ(θ, x), as is the case for unpolarised Compton scattering. The right panel of Fig. <ref> visualises the scattering phase function for three (incoming) polarisation states, illustrating how the phase function for unpolarised photons (P_L=0) is symmetric around the incoming photon direction , while for fully polarised photons (P_L=1), the phase function is almost perfectly symmetric around the incoming polarisation direction . For partially polarised photons, the scattering phase function is truly asymmetric, as illustrated for P_L=0.5 on the right panel of Fig. <ref>.

After a Compton scattering interaction, the Stokes parameters of the scattered photon are updated as

    [ I'; Q'; U'; V' ]∝[ S_11(θ, x)   I + S_12(θ, x)  (cos2φ Q+sin2φ U);  S_12(θ, x)   I + S_22(θ, x) (cos2φ Q+sin2φ U);                  S_33(θ, x) (-sin2φ Q+cos2φ U);                                   S_44(θ, x) V ],

combining Eq. (<ref>), (<ref>), and (<ref>). In addition, the photon energy is reduced to E'=C(θ, x) × E, and the updated Stokes vector ' is normalised so that I'=C(θ, x) × I, to comply with the conservation of four-momentum for the photon-electron pair.


Through Compton scattering, the radiation field obtains a linear polarisation degree P'_L, which depends on the initial polarisation state (P_L and γ) and the scattering direction (θ and φ):

    P'_L = √((  S_12  + S_22 P_Lcos2(γ-φ))^2 + (S_33 P_Lsin2(γ-φ))^2)/S_11 + S_12P_Lcos2(γ-φ),

combing Eq. (<ref>), (<ref>) and (<ref>), with S_ij≡S_ij(θ, x) being the Müller matrix elements given by Eq. (<ref>). For various incident polarisation degrees P_L, the resulting polarisation degree P'_L is shown in Fig. <ref>, illustrating how any polarisation degree between 0 and 1 can be obtained, depending on the specific scattering direction. Regardless of the initial polarisation state, scattered photons are maximally polarised (P'_L=1) for θ=90^∘, while forward scattering and backscattering leaves the polarisation degree unchanged. The resulting polarisation degree is shown at 2 keV (solid line) and 8 keV (dashed line), demonstrating how the energy dependence of P'_L is negligible in the IXPE range. 

Equivalently, the linear polarisation angle obtained through Compton scattering is

    γ' = 1/2arctan_2(S_33(θ, x)  P_Lsin2(γ-φ)/S_12(θ, x) + S_22(θ, x)  P_Lcos2(γ-φ)),

combing Eq. (<ref>), (<ref>) and (<ref>). However, we remind that γ' as given by Eq. (<ref>) refers to a different reference direction for each pair of scattering angles (θ, φ) (see Eq. (<ref>)). Therefore, these polarisation angles cannot be directly compared, except for some special symmetric cases such as the scenario discussed in Sect. <ref>. For unpolarised photons (P_L=0), we recover the standard result that scattering induces polarisation perpendicular to the scattering plane (γ'=90^∘, with ' in the scattering plane). 



§ SKIRT IMPLEMENTATION




 §.§ The radiative transfer code SKIRT


For the remainder of this work, we focus on radiative transfer simulations with the radiative transfer code SKIRT <cit.>. SKIRT is a state-of-the-art Monte Carlo radiative transfer code, developed and maintained at Ghent University, which implements a Monte Carlo photon life cycle emulating absorption, scattering, and re-emission in complex 3D transfer media <cit.>. SKIRT simulations in full 3D are facilitated by the implementation of various acceleration techniques <cit.>, advanced grids for discretising the transfer medium <cit.>, and a hybrid parallelisation scheme which combines multi-threading with multi-processing <cit.>. The SKIRT code is open-source, well-documented[<https://skirt.ugent.be>], and publicly available online[<https://github.com/SKIRT/SKIRT9>], with tutorials for users and developers.

SKIRT models absorption, scattering, and re-emission in dusty astrophysical systems <cit.>, which includes emission from stochastically heated dust grains <cit.>, polarisation due to scattering on spherical dust grains <cit.>, and polarised emission from aligned spheroidal dust grains <cit.>. Beyond dust, SKIRT models scattering on free electrons <cit.>, absorption and emission at the 21 cm line of neutral hydrogen <cit.>, non-LTE line radiative transfer in the (sub)mm and infrared <cit.>, and resonant line scattering of H Lyα photons <cit.>. Recently, the SKIRT code was extended into the X-ray regime, to model Compton scattering on free electrons, photo-absorption and fluorescence by cold atomic gas, scattering on bound electrons, and extinction by dust <cit.>. Furthermore, the kinematics of moving sources and moving transfer media are self-consistently incorporated into the SKIRT radiative transfer calculations, so that the effect of bulk velocities and velocity dispersions can be properly modelled <cit.>.

The SKIRT code features a large suite of  model geometries, radiation sources, medium characterisations, instruments, and probes <cit.>, in addition to interfaces for post-processing hydrodynamical simulations <cit.>. SKIRT can calculate self-consistent fluxes, images, spectra and polarisation maps from mm to X-ray wavelengths, with recent applications in  galaxies <cit.>, active galactic nuclei <cit.>, and other fields of astrophysics <cit.>.



 §.§ X-ray polarisation implementation


The polarisation framework for tracing Stokes vectors in SKIRT has been implemented by <cit.>, while the X-ray physical processes in cold neutral media (without polarisation support) have been implemented by <cit.>. As photo-absorption and fluorescence do not depend on the polarisation state of the incoming photon, the original implementation is retained for these processes, with the important distinction that the Stokes vector is reset after fluorescent re-emission (i.e. the Q', U', and V' parameters are set to zero; see Sect. <ref>).

In this work, bound-electron scattering in cold atomic gas is approximated as free-electron scattering on Z free electrons per neutral gas atom, as discussed in Sect. <ref>. For a gas column with column density N_H, the optical depth for scattering is thus

    τ(E) = (∑_Z  a_Z  Z) · N_H·σ(E) ,

with a_Z being the number density of element Z relative to H i, and σ(E) the Compton scattering cross-section given by Eq. (<ref>).
As this optical depth does not depend on the polarisation state of the incoming photon, sampling for scattering positions is conceptually identical to the free-electron scattering implementation as presented in <cit.>. Given the abundance table { a_Z}, one can easily generate a random interaction point in the forward  direction, following for example <cit.>. 

Once the scattering location is set, a pair of scattering angles (θ, φ) can be sampled from the phase function Eq. (<ref>) using the conditional distribution method. Due to its symmetry around Ψ(θ, x), the azimuthal modulation term does not contribute to the marginal probability distribution for θ, and therefore, random scattering angles θ can be generated from the same univariate distribution p(θ; x) as for unpolarised Compton scattering:

    p(θ; x) 
           = 3 σ_T S_11(θ, x)sinθ/8 σ(x).

This distribution for θ is a complex function of the incoming photon energy, but can be sampled efficiently using a variation on Khan's technique described by <cit.>. Once a scattering angle θ is generated, we obtain the conditional distribution for the azimuthal scattering angle φ as

    p_θ(φ; x, S) 
       = 1/2π( 1 + S_12(θ, x)/S_11(θ, x)P_Lcos2(φ-γ)),

which explicitly depends on the incoming polarisation state. A random scattering angle φ can be sampled from Eq. (<ref>) through numerical inversion as

    χ   = ∫_0^φ p_θ(φ'; x, S)  dφ'
       = 1/2π( φ + S_12(θ, x)/S_11(θ, x)P_Lsinφcos(φ-2γ)),

with χ being a uniform deviate between 0 and 1.

Once the scattering angles θ and φ are generated, the photon energy of the interacting photon is updated to E'=C(θ, x) × E and the Stokes parameters are updated as described by Eq. (<ref>) (normalised so that I'=C(θ, x) × I, as discussed in Sect. <ref>). Hereafter, the scattered photon continues its way through the transfer medium, in a direction ' that can be calculated from the original direction  and the scattering angles (θ, φ).

Eventually, all photons packets are recorded by (a set of) SKIRT instruments (corresponding to preset observer directions), using the peel-off technique as described in <cit.>. We remind that the Stokes vectors of recorded photons are first rotated to correspond to the observer's north direction (following the IAU conventions; see Sect. <ref>), before they are binned to produce Stokes spectra and polarisation maps.



 §.§ Implementation verification




  §.§.§ Simulation setup



To verify the X-ray polarisation implementation in SKIRT, we set up dedicated SKIRT simulations, specifically designed to recover the details of individual photon-electron interactions. The simulation setup is shown in Fig. <ref> and consists of a collimated beam illuminating a cloud of free electrons at the origin. The 4 keV source photons are emitted in the positive z-direction  and are linearly polarised along the x-axis, with P_L=0.5. Without loss of generality, we can choose =e_x as the reference direction (with e_x being the Cartesian unit vector in the positive x-direction), so that γ = 0 by definition. The central electron cloud is small compared to its distance to the source and has a column density of 10^19 cm^-2, making the flux contribution of photons that scatter more than once negligible (i.e. <0.001% of the one-time-scattered flux). The distance to the observer is D, and the spherical coordinates of the observer direction ' are just the scattering angle θ and the azimuthal scattering angle φ.

For any observer direction ', the scattered photon flux can be obtained with SKIRT, and because of the peel-off detection technique, this is possible for any exact (θ, φ)-direction, with no spurious blurring due to averaging over some finite solid angle around ' <cit.>. Furthermore, the `smart' photon detectors in SKIRT can recover the individual flux contributions of direct (i.e. non-interacting) and reprocessed photons <cit.>, in addition to recording the total Stokes I, Q, U, and V spectra. For all (non-forward) directions, the total observed flux is equal to the scattered flux component, which is virtually equal to the one-time-scattered flux  (as further interactions are negligible at N_e = 10^19 cm^-2).

As a first sanity check, we calculated the ratio of the observed direct spectrum over the input spectrum for the forward direction (θ=0), recovering the expected result of 1-exp[-N_e·σ(E)], with σ(E) being the scattering cross-section given by Eq. (<ref>). In the next subsections, we use the SKIRT simulation output to recover the phase function (Sect. <ref>), the polarisation degree (Sect. <ref>), and the polarisation angle (Sect. <ref>), which are then verified against the analytical formula obtained in Sect. <ref>.



  §.§.§ Phase function


The scattering phase function Φ for polarised photons (P_L>0) depends explicitly on the azimuthal scattering angle φ, which is not the case for unpolarised Compton scattering (see Eq. (<ref>)). Using the SKIRT output for the simulation setup described in Sect. <ref>, we can infer the scattering phase function Φ_SKIRT as implemented in SKIRT by relating the (one-time) scattered photon flux density I(θ, φ) in a direction ' to the fraction of the specific beam luminosity L that interacted inside the electron cloud:

    Φ_SKIRT(θ, φ; E, )= I(θ, φ; E')· C(θ, E)^2 · D^2/(1-exp[-N_e·σ(E)])· L(E),

with I(θ, φ) being the photon flux measured at E' = C(θ, E) × E, as Compton scattering is inelastic, and L and D being known input parameters. We note that the factor C(θ, E)^2 appears in the nominator of Eq. (<ref>) as the interval dE around E corresponds to the interval dE'=C(θ, E)^2×dE around E' after scattering.


Fig. <ref> shows the scattering phase function Φ_SKIRT, which was inferred from the SKIRT results for the simulation setup described in Sect. <ref>. Comparing this Φ_SKIRT to the theoretical phase function given by Eq. (<ref>), we find an excellent agreement for all scattering angles (θ, φ), recovering the prominent azimuthal modulation that governs the scattering physics for polarised photons. Using Eq. (<ref>), we verified the SKIRT results for a range of polarisation states  and photon energies E (especially beyond the IXPE range, where the energy dependence becomes more pronounced), assuring a correct implementation of polarisation-dependent Compton scattering in SKIRT.



  §.§.§ Polarisation degree



The Stokes I(θ, φ), Q(θ, φ), and U(θ, φ) spectra calculated with SKIRT can be used to obtain the polarisation degree for the setup described in Sect. <ref> by applying Eq. (<ref>). Fig. <ref> shows the inferred polarisation degree as a function of the scattering angle, together with the theoretical formula for P_L', given by Eq. (<ref>). As the polarisation degree induced through Compton scattering is virtually symmetric around θ=90^∘ (see Fig. <ref>), we focus on scattering angles θ≤ 90^∘ only. We find an excellent agreement between the SKIRT results and the theoretical formula Eq. (<ref>), recovering the complex behaviour of P_L' with θ and φ. We find that forward scattering (θ= 0) does not change the polarisation degree, while P_L'=1 when θ=90^∘, as predicted by Fig. <ref>.



  §.§.§ Polarisation angle



Similarly, the linear polarisation angle γ' can be calculated from the Stokes spectra obtained with SKIRT by applying Eq. (<ref>). The inferred polarisation angle γ'_SKIRT is then defined relative to the north direction of the observer frame, as discussed in Sect. <ref>. This γ'_SKIRT can be directly compared to the theoretical formula as given by Eq. (<ref>), as for this particular setup (Sect. <ref>), Eq. (<ref>) refers to a reference direction ' that is also the observer north direction[Using Eq. (<ref>) with =e_z and =e_x (see Sect. <ref>), ' is just equal to the spherical unit vector e_θ. The observer north direction is -e_θ, but Stokes vectors are invariant under rotations over 180^∘ (see Eq. (<ref>)).].

The linear polarisation angle inferred from the SKIRT output is shown in Fig. <ref>, for various scattering angles. For forward scattering and backscattering (θ= 0 and 180^∘, respectively), the observed γ' is just the projected[Projected on the sky as seen from the observer location.] incident polarisation direction, as the polarisation direction is retained by scattering at θ = 0^∘ or θ = 180^∘ (see Eq. (<ref>) and (<ref>)). For intermediate scattering angles, the observed polarisation angle lays between γ'=90^∘ (i.e. perpendicular to the scattering plane) and the specific angle corresponding to the projected incoming polarisation direction. For φ = 90^∘ or 270^∘, the projected incident polarisation direction is zero, and therefore, the observed polarisation direction is just perpendicular to the scattering plane (γ'=90^∘),  identical to the case of unpolarised incident photons (see Sect. <ref>). 

We find an excellent agreement between the SKIRT results and the theoretical formula Eq. (<ref>) for γ', recovering the complex behaviour of the polarisation angle with the scattering direction. Together with the results of Sect. <ref> and Sect. <ref>, this assures us that the Compton scattering formula are implemented correctly in SKIRT.



§ TORUS MODELS




 §.§ Model setup


As a first application of the new X-ray polarisation capabilities of the SKIRT code, we explored the spectro-polarimetric properties of a 2D toroidal reprocessor of cold gas, modelling the parsec-scale circumnuclear medium of AGN. This medium is expected to reprocess the primary X-ray emission of the central X-ray corona, producing a distinct polarisation signal which can be used to constrain the geometry of the reprocessor in heavily obscured AGN <cit.>. 


We adopted a wedge torus geometry of uniform-density gas, centred on an isotropic point source representing the X-ray corona (see Fig. <ref>). This torus geometry is identical to the torus geometries of the BNTorus <cit.> and borus <cit.> spectral models and is also known as a `flared disc'. The X-ray radiative transfer problem is scale invariant, and therefore the model geometry is fully defined by the torus covering factor CF = cosθ and the equatorial hydrogen column density N_H, which are indicated on Fig. <ref>. For the wedge torus geometry, all obscured sightlines (cosi < CF) have the same line-of-sight N_H, while unobscured sightlines have N_H = 0.

We considered photo-absorption, fluorescence, and scattering by cold neutral gas as described in Sect. <ref>, and  adopted <cit.> solar gas abundances. We assumed a standard power law spectrum for the central X-ray source, characterised by a photon index Γ and an exponential cut-off energy E_cut. The primary X-ray photons were assumed to be unpolarised (P_L=0), so that the emerging polarisation signal could be entirely attributed to reprocessing in the AGN torus. Indeed, IXPE observations of unobscured AGN indicate that the coronal polarisation levels are low <cit.>, so that the initial polarisation state of the coronal X-ray photons is easily washed out by multiple scattering in the AGN torus <cit.>.

We set up radiative transfer simulations in the wedge torus geometry, using the SKIRT code with X-ray polarisation enabled[SKIRT version 9 git commit .]. First, we specified the spatial grid on which the torus medium was discretised. In this case, the axial symmetry of the torus could be exploited to build a 2D grid, which would speed up the radiative transfer simulations significantly. Furthermore, the density distribution in the azimuthal plane could be gridded using a 2D polar grid with angular bins coinciding with the torus edges, so that the gridded distribution would equal the exact model distribution, with no discretisation effects.

Stokes spectra were calculated over the 0.3 to 200 keV range, adopting 1000 logarithmic wavelength bins, which corresponds to a spectral resolution of R=154. In SKIRT, X-ray radiative transfer interactions were modelled up to 500 keV, as the X-ray reflection spectrum and the corresponding polarisation signal are produced by Compton down-scattering at energies beyond 200 keV <cit.>.



 §.§ Radiative transfer results



The SKIRT radiative transfer results are shown in Fig. <ref> for one particular realisation of the wedge torus model (CF=0.45 and cosi=0.4), with different panels representing different values for the column density, ranging from logN_H=23.4 (Compton-thin, top left) to logN_H=25.0 (Compton-thick, bottom right). The adopted source parameters are representative for the nearby AGN in the Circinus galaxy (Γ = 1.8,  E_cut=200 keV, and L_2-10=2.8 × 10^42 erg s^-1, observed at a distance of 4.2 Mpc). We focus on the SKIRT results in the 2-8 keV IXPE range.

The Stokes I, Q, and U spectra are shown at the top of each panel. The observed Stokes U parameter is always zero, which was expected as the north direction of the observer frame was chosen to just coincide with the projected symmetry axis of the torus (while further considering that a 2D system cannot exhibit features that are asymmetric relative to its projected symmetry axis). In a realistic observational context, the observed system will be rotated about the line of sight of the observer, mixing the Q and U parameters as described by Eq. (<ref>). This rotation-induced mixing of the Stokes parameters can be observed, which could be used to infer the orientation of the principal axis of the AGN system, forming a powerful probe on the (spatially unresolved) circumnuclear medium <cit.>.

From these Stokes spectra, the linear polarisation degree P_L and the linear polarisation angle γ were calculated as a function of the photon energy, using Eq. (<ref>) and (<ref>). As Stokes U is always zero, γ can only be 0^∘ or 90^∘. Therefore, one could visualise the observed polarisation direction with positive P_L values denoting polarisation in the north direction (γ=0^∘) and negative P_L values denoting polarisation in the horizontal direction (γ=90^∘) (see P_L, in purple, at the bottom of each panel in Fig. <ref>). As the north direction of the observer frame was chosen to just coincide with the projected torus axis, γ=0^∘ and γ=90^∘ correspond to polarisation parallel and perpendicular to the torus symmetry axis, respectively.

The radiative transfer results shown in Fig. <ref> demonstrate the complex behaviour of the polarisation observables of a 2D torus model in the 2-8 keV range. Most prominently, we note that the Stokes Q spectra have a fixed negative sign up to a certain photon energy E_flip (indicated in green on Fig. <ref>), beyond which Q becomes positive. With Stokes U being zero, this means that for E<E_flip, the observed emission is polarised perpendicular to the torus symmetry axis (i.e. γ=90^∘), while at higher energies, the polarisation is parallel to this symmetry axis (i.e. γ=0^∘).

The exact photon energy E_flip where the polarisation angle γ flips from 90^∘ to 0^∘ is a function of the opacity of the torus medium, increasing from 2.0 keV at logN_H=23.4 to 6.6 keV at logN_H=24.8[For logN_H=25.0, E_flip = 10.4 keV, outside of the IXPE range.]. Furthermore, the effect of the torus opacity is also observed in the logN_H=24.8 panel of Fig. <ref>, where the polarisation angle γ flips back to 90^∘ at 7.1 keV  (i.e. at the Fe K absorption edge), where the torus opacity rises discontinuously. A second flip from γ=90^∘ to γ=0^∘ is then observed at 8.8 keV. This behaviour of the polarisation angle with photon energy is discussed in Sect. <ref>.

Away from E_flip, Stokes Q is relatively featureless and roughly follows the spectral shape of the reflected continuum. In particular, the  polarised flux (i.e. Stokes Q) does not contain any fluorescent lines, as fluorescent line photons are emitted at P_L=0 (see Sect. <ref>). However, as the total flux (i.e. Stokes I) contains strong fluorescent lines, the linear polarisation degree is heavily diluted at the fluorescent line energies, which produces line features where P_L approaches zero (see Fig. <ref>). Similarly, at low N_H, the polarisation signal is heavily diluted by the direct flux of the (unpolarised) primary X-ray source.



 §.§ Model grid



We calculated a grid of AGN torus models for observational data fitting based on the wedge torus geometry described in Sect. <ref> by varying N_H between 10^22 and 10^26 cm^-2, the photon index Γ between 1 and 3, and the exponential cut-off energy log E_cut between 1.5 and 2.9 (corresponding to 30 keV and 800 keV, respectively). We considered torus covering factors between 0.05 and 0.95, observed at inclinations between 0^∘ and 90^∘, with all parameters being sampled as described in Table <ref>. This leads to 203280 unique torus model realisations, forming a parameter space that is covering most of the obscured AGN observed in the local Universe <cit.>. For each parameter combination, the Stokes I, Q, and U spectra were calculated from 0.3 to 200 keV, with a spectral resolution of R=154 (see Sect. <ref>).

For all model realisations, the number of simulated photon packets was kept at 10^8, resulting in output spectra with limited Monte Carlo noise. The simulation run time of a single model realisation depends mostly on the covering factor and the column density of the torus, as more material requires more interactions to be modelled. In addition, the spectral hardness of the source has a small effect, as hard X-ray photons experience more scattering events before being absorbed. For a standard source spectrum with Γ = 1.8 and E_cut=200 keV, individual simulations take between 0.9 and 33.1 minutes on a modest 2.2 GHz 16-core node, demonstrating the computational efficiency of the SKIRT code. By using the high-performance computing infrastructure at the Flemish Supercomputer Centre, it was possible to run all 203280 torus model realisations in just 55 hours.

The radiative transfer results were converted to an XSPEC <cit.> table model named xskirtorsmooth, which is publicly released with this work[<https://github.com/BertVdM/xskirtor>]. This model can be used for observational data fitting within XSPEC as

    model polrot × atable{xskirtor smooth.mod}.

In addition to the five free model parameters listed in Table <ref>, the table model requires a redshift and a luminosity normalisation[Following the same convention as the cutoffpl model in XSPEC, the norm parameter is defined as the constant factor K in the unobscured flux density of the primary source: F_E(E) = K  E^-Γexp(-E/E_cut), in units of counts/s/keV/cm^2. The norm parameter is approximately equal to the unobscured flux density at 1 keV.], providing the physical scaling of the model spectra. polrot then sets the roll angle θ of the torus system around the line of sight of the observer, as the table model corresponds to an observer that has its north direction coinciding with the projected torus symmetry axis. This brings the number of free model parameters to eight. The xskirtorsmooth model can be further combined with other XSPEC models, such as a tbabs component to model galactic foreground extinction <cit.>.

The xskirtorsmooth model predicts both X-ray spectra and X-ray spectro-polarimetry, enabling the simultaneous fitting of spectral and polarimetric observations in the X-ray band. As both the spectral coverage and the spectral resolution of the model templates are largely exceeding the capabilities of modern X-ray polarisation observatories, they can be applied for the interpretation of observational IXPE, XPoSat, and eXTP data, proposal writing, and the definition of future missions.

The corresponding spectral model (modelling Stokes I) is well suited for fitting CCD-based X-ray spectra as obtained with modern X-ray observatories such as XMM-Newton, Chandra, NuSTAR, Swift, INTEGRAL, AstroSat, and more. Furthermore, we provided an additional XSPEC table model for the 1.5 to 15 keV subrange, with an adaptive spectral resolution achieving Δ E = 0.5 eV around the strongest fluorescent lines, for fitting the high-resolution X-ray spectra that are being obtained with the XRISM/Resolve <cit.> and will be obtained with the Athena/X-IFU <cit.>.



§ DISCUSSION




 §.§ Polarisation angle as a function of the photon energy



In Sect. <ref>, we described how the reprocessed torus emission is polarised perpendicular to the torus symmetry axis (γ=90^∘) at low energies, while it is parallel to the torus symmetry axis (γ=0^∘) at higher energies. This behaviour of the polarisation angle γ with photon energy is observed for a wide range of torus parameter combinations (see Sect. <ref>), with a transition energy E_flip that scales with the torus opacity. To understand how this behaviour is related to the torus geometry, we can inspect the spatially resolved Stokes surface brightness maps, which can be calculated with SKIRT for each parameter combination.

As a demonstration, we focus on the logN_H=23.8 torus model realisation that was discussed in Sect. <ref>, for which the Stokes spectra are shown in the top right panel of Fig. <ref>. The corresponding Stokes Q (left), Stokes U (middle), and Stokes I (right) surface brightness maps are shown in Fig <ref>, with the top row representing the 2 keV to E_flip energy band, and the bottom row representing the E_flip to 8 keV band, with E_flip = 2.9 keV. In addition, the total surface brightness map (Stokes I, on the right) is overlaid with a linear polarisation map, visualising the polarisation degree and polarisation angle as derived from the Stokes Q and Stokes U surface brightness maps.

Looking at the top right panel of Fig. <ref>, we find that for E<E_flip, the total observed flux is dominated by reprocessed photons that scattered off the illuminated (unobscured) backside of the torus. For E>E_flip on the other hand (bottom right), source photons are able to penetrate the obscuring front side of the torus, so that the total observed flux is dominated by direct source emission. However, this unpolarised direct flux merely dilutes the observed polarisation signal, without influencing the polarisation angle γ that is eventually observed. Therefore, we can focus on the reprocessed flux only, finding that for E>E_flip, scattered photons originate from a more extended region on the illuminated backside of the torus, while they can also escape through the obscuring front side of the torus, as evident from the bottom right panel of Fig. <ref>.

Using the smart photon detectors in SKIRT (see Sect. <ref>), we can confirm that the reprocessed torus emission is dominated by one-time-scattered photons, which significantly simplifies the interpretation of the linear polarisation maps shown in Fig. <ref>. As the primary source photons are unpolarised (see Sect. <ref>), Compton scattering induces polarisation that is exactly perpendicular to the scattering plane (see Eq. (<ref>)). This means that the projected polarisation pattern is always circular (i.e. in each pixel, the observed polarisation direction is perpendicular to the direction towards the central pixel), as shown in Fig. <ref>. Indeed, the Stokes Q surface brightness is positive in the east and west quadrants (describing polarisation in the vertical direction) and negative in the north and south quadrants (describing polarisation in the horizontal direction). Similarly, Stokes U is positive in the north-east and south-west quadrants and negative in the north-west and south-east quadrants (see Fig. <ref>). The polarisation direction of individual reprocessed photons thus depends on the specific sky region of the final scattering interaction inside the torus (in this specific case, dominated by single scattering). 

In addition, we find that also the observed polarisation degree depends on the sky region, with P_L being high in the east and west quadrants and low in the north and south quadrants (see Fig. <ref>, where P_L is visualised by the length of the polarisation map segments). This behaviour is a direct result of the projected inclined torus geometry: In the north and south quadrants, the distribution of last scattering angles (i.e. to reach the observer) forms a narrow peak centred on 130^∘ and 50^∘, respectively. For both quadrants, this correspond to a distribution of polarisation degrees having P_L<0.3 for 60% of all photons and P_L<0.5 for 90% of all photons, as predicted by Eq. (<ref>). On the other hand, the distribution of scattering angles is much broader in the east and west quadrants and peaks at 90^∘, which corresponds to higher average polarisation degrees (see Fig. <ref>). Indeed, in these side quadrants, 50% of all photons are having P_L>0.6.

Summarising, the net (i.e. spatially integrated) polarisation direction that is eventually observed, is the result of the precise balance between the polarised flux originating from the different sky regions, with both the total flux and the polarisation degree of each sky region being closely related to the torus geometry. As the projected torus geometry is perfectly symmetric around the observer north direction, the Stokes U surface brightness maps are perfectly antisymmetric around this axis (see Fig. <ref>), so that the spatially integrated Stokes U fluxes are always zero, as discussed in Sect. <ref>. The Stokes Q surface brightness maps on the other hand do not exhibit such trivial symmetry, so that the spatially integrated polarisation angle γ depends on the balance between the polarised flux in the east and west quadrants (where Stokes Q>0 and γ=0^∘) compared to the north and south quadrants (where Stokes Q<0 and γ=90^∘).

We can now explain the behaviour of the polarisation angle with photon energy for this particular setup: At low energies (E<E_flip), the reprocessed torus emission is purely dominated by photons that scattered off a small region on the illuminated backside of the torus. This backside region mostly covers the northern sky quadrant where Stokes Q is negative, so that the spatially integrated polarisation signal is perpendicular to the torus symmetry axis (γ=90^∘). At higher energies (E>E_flip), the reprocessed torus emission originates from a more extended region on the torus backside and also escapes through the obscuring front side of the torus, so that all four sky quadrants are covered. At E_flip, the polarised flux is still dominated by the unobscured torus backside, but now the positive Q contributions on the torus backside (close to the torus front edge, as shown in red) start to dominate over the negative Q contributions (shown in blue; see the bottom left panel of Fig. <ref>). Indeed, while most of the reprocessed flux is still contained within the northern sky quadrant, we find that the polarisation degree is significantly higher in the east and west quadrants, so that the polarised flux is actually dominated by the contribution of these side quadrants, and the net polarisation is parallel to the torus axis (γ=0^∘). At even higher energies, the polarised flux is dominated by photons escaping thought the obscuring front side of the torus, mostly from the east and west sky regions, so that Q>0 and γ=0^∘.

The flux balance between the northern sky region (Q<0) and the east and west sky regions (Q>0) is determined by the level of obscuration of the latter regions. This explains why the transition energy E_flip scales with the torus opacity (see Fig. <ref>), as higher torus column densities require higher photon energies (with more penetrating power) to escape from the east and west sky regions. This is a direct result of the specific torus geometry, and we conclude that spatially resolved Stokes surface brightness maps form a powerful tool to study the geometrical effects that are encoded in spectro-polarimetric observations. The SKIRT code allows for calculating these surface brightness maps at a high signal-to-noise, in limited computational time.



 §.§ Inclination - covering factor contours



Inspired by Sect. <ref>, where we found an interesting behaviour of the polarisation observables at cos i ≲CF, this section discusses the wedge torus results of Sect. <ref> as a function of the cos i and CF model parameters. This analysis should also be useful for the interpretation of observational data, to constrain torus properties from broadband polarimetry. Fig. <ref> shows the total polarisation degree (over the 2-8 keV band) as a function of cos i and CF, with different panels representing different values for the torus column density (increasing from logN_H=23 to logN_H=25).

We find that the total polarisation degree is virtually zero (P_L<0.3%) for unobscured sightlines (cos i > CF), where the polarisation signal of the AGN torus is mostly diluted by the direct source emission[This is when assuming unpolarised primary emission (Sect. <ref>). Alternatively, the torus signal is `diluted' by the polarisation signal of the X-ray corona <cit.>.]. The polarisation degree reaches a maximum at cos i ≲CF, where the illuminated backside of the torus can be observed without significant obscuration, similar to Sect. <ref>. Furthermore, the total polarisation degree is observed to increase from 1% to 30% when logN_H increases from 23 to 25, as more torus material results in more scattering interactions inducing a stronger polarisation signal, while also the unpolarised direct flux component is more obscured.


Finally, two distinct polarisation maxima can be observed in each panel of Fig. <ref>, separated by a region where P_L approaches zero. One local maximum is located at cos i ≲CF = 0.7  to  0.9 (top right corner of each panel, with γ=90^∘), while the other local maximum is located at cos i ≲CF = 0.2  to  0.4 (bottom left corner of each panel, with γ=0^∘). Similar to Sect. <ref>, we calculated the spatially resolved Stokes surface maps to study this behaviour, which are shown in Fig. <ref>. We find that both of these polarisation maxima are related to scattering on the illuminated backside of the torus, similar to Sect. <ref>.

At cos i ≲CF = 0.7  to  0.9 (i.e. tori that are almost entirely closed), the illuminated backside of the torus mostly covers the northern sky quadrant (where Stokes Q<0), so that the net polarisation signal is perpendicular to the torus axis (γ=90^∘). For cos i ≲CF = 0.2  to  0.4 on the other hand (i.e. thin, disc-like tori), the illuminated backside mostly covers the east and west quadrants (where Stokes Q>0), so that the polarisation is parallel to the torus axis (γ=0^∘). For parameter combinations in between these two maxima, the torus backside covers all three sky quadrant, so that positive and negative Stokes Q contributions (partially) cancel out, eventually reaching P_L=0 at the border region that separates the two local maxima.



 §.§ Polarisation angle flip for different sightlines


In Sect. <ref>, we presented Stokes spectra for one sightline only: cos i ≲CF, which offers an unobscured view on the illuminated backside of the torus. The reason for this specific choice is that we plan to compare these results to simulations that include a polar component in follow-up work, showing similar spectra[Indeed, polar-extended dusty gas could act as an illuminated reflector that is visible through unobscured sightlines, which would require less fine-tuning of the observer inclination (forming a natural solution for strong reflection spectra).]. However, because of this very specific viewing angle, the trends that were found for CF=0.45 and cos i =0.4 in Sect. <ref> might differ from more general trends at arbitrary obscured sightlines, which we investigate in this section. In particular, we focus on the photon energy E_flip where the polarisation angle γ flips from 90^∘ to 0^∘, as a function of the torus column density.


Using the radiative transfer results of the torus model grid that was presented in Sect. <ref>, we measured E_flip for a broad range of torus parameter combinations, focussing on the same source parameters as Sect. <ref> (Γ = 1.8 and E_cut=200 keV), but a higher torus covering factor (CF=0.85) to allow for more obscured sightlines (cos i < CF). In Fig. <ref>, the polarisation flip energy E_flip is shown as a function of the torus column density for different (obscured) observer inclinations, generalising the trend that was found in Sect. <ref>: A linear relation between logN_H and logE_flip is found for all viewing angles, which is broken at the Fe K absorption edge at 7.1 keV. Indeed, similar to Sect. <ref>, E_flip is found to correlate with the torus opacity, which increases discontinuously at the Fe K edge for a fixed column density. 

For CF=0.85, the polarisation flip happens when the torus medium becomes sufficiently transparent, so that the reprocessed flux escaping through the front side of the torus (mostly Q>0) starts to dominate over the flux related to scattering on the torus backside (Q<0). As the former region is obscured by the torus while the latter region is not, the ratio between these two components naturally scales with the torus opacity, explaining the trend with logN_H (see Fig. <ref>). In addition, as the flux of the torus backside decreases with increasing inclination (which is a projection effect), the transmitted flux (having Q>0) starts to dominate at lower E_flip for higher inclinations (i.e. lower cos i), explaining the observed trend with cos i in Fig. <ref>.

Finally, we note how the details of the polarisation angle flip depend on the specific combination of CF and cos i. For CF=0.85, we found a balance between transmission through the torus front side and scattering on the torus backside for all cos i, which is different from the scenario discussed in Sect. <ref>, where we found a balance between positive and negative Stokes Q contributions on the torus backside for CF=0.45. For smaller covering factors, the positive Stokes Q contributions on the torus backside could be dominating over the entire IXPE range (so that no polarisation flip is observed), while for some CF and cos i combinations, two flips could occur. These intricacies, and their link to the torus geometry, can be studied in great detail with SKIRT, in particular based on Stokes surface brightness maps, which can be calculated in little computational time.



 §.§ XSPEC torus model


With this work, we release an AGN torus model that describes both X-ray spectra and X-ray polarisation for observational data fitting with XSPEC (see Sect. <ref>). This xskirtorsmooth model represents a smooth toroidal reprocessor of cold gas, positioned in the equatorial plane, as described in Sect. <ref> (see Fig. <ref>). While similar smooth torus models have been very successful in modelling the observational X-ray spectra of obscured AGN <cit.>, these models do not incorporate geometrical complexities such as clumpy or filamentary substructures, or polar-extended dusty gas, which might be omnipresent in local AGN (Sect. <ref>). Therefore, the xskirtorsmooth model should be used as a tool to gain insights into the representative properties of the circumnuclear medium (such as its covering factor or average column density), more than the final solution for the AGN torus geometry, when being applied to observational data. 

The xskirtorsmooth model is provided in two different flavour variations: a coupled configuration which provides the direct and reprocessed flux components as a single table, and a decoupled configuration which allows for varying the line-of-sight column density independently from the equatorial column density[In fact, all model parameters could be varied independently (given that there would be a physical motivation to do so).] <cit.>. As a final remark, we note that the exponential cut-off energy of the primary X-ray source (Sect. <ref>) has a noticeable effect in the 2-8 keV range, even when E_cut>100 keV[For example, for E_cut = 100 keV, the effect of the exponential cutoff on the X-ray spectrum at 8 keV is 8%.]. Therefore, the logE_cut model parameter (see Table <ref>) is relevant beyond the modelling of hard X-ray spectra and should not be neglected at lower X-ray energies such as the IXPE range.



§ SUMMARY AND OUTLOOK


In this work, we presented a general framework for modelling X-ray polarisation in 3D radiative transfer simulations of cold gas and dust, to model the spectro-polarimetric signal that is produced by X-ray reprocessing in AGN circumnuclear media. We discussed how radiative transfer processes depend on the polarisation state of the incoming photon and how the polarisation state is updated by Compton scattering and fluorescence (Sect. <ref>). We described how polarised X-ray radiative transfer can be implemented using a Monte Carlo method and provided an implementation to the 3D radiative transfer code SKIRT, which is publicly available online[<https://skirt.ugent.be>] (Sect. <ref>).

As a first application, we focussed on a 2D torus geometry in Sect. <ref>, to demonstrate the new X-ray polarisation capabilities of the SKIRT code without going into the details of 3D structure and its effect on spectro-polarimetric observables. However, we note that the current SKIRT implementation works in full 3D already, so that more complex models <cit.> can be run with X-ray polarisation enabled (i.e. without any further modifications to the code). In future work, we will focus on these models with a truly 3D structure beyond the classical torus.

For the 2D wedge torus model (Sect. <ref>), we calculated Stokes spectra at a high signal-to-noise (Fig. <ref>) and computed the linear polarisation angle and polarisation degree as a function of photon energy. We found that the polarisation angle flips from 90^∘ to 0^∘ at a specific energy inside the IXPE range (Sect. <ref>), which we interpreted as a balance between the reprocessed flux originating from different regions of the torus, with a direct link to the torus geometry (Sect. <ref>). Furthermore, we found that the polarisation degree reaches a maximum at (obscured) sightlines having cos i ≲CF (Sect. <ref>), where the torus backside can be observed without significant obscuration. However, depending on the torus covering factor, this polarisation maximum can be parallel or perpendicular to the torus axis (see Fig. <ref>). Finally, we found that the specific photon energy E_flip where the polarisation angle γ flips from 90^∘ to 0^∘ scales with the torus opacity (Sect. <ref>), as the polarisation flip is related to a balance between torus regions with different levels of obscuration. These intricacies, and their link to the torus geometry, were studied in great detail using the Stokes surface brightness maps, which SKIRT can calculate in a short amount of computational time.

With this work, we release spectro-polarimetric templates for fitting observational data of obscured AGN based on the torus model grid presented in Sect. <ref> (and discussed in Sect. <ref>). This  X-ray torus model is provided as an XSPEC table named xskirtorsmooth, which can simultaneously describe X-ray spectra and spectro-polarimetry over the 0.3 to 200 keV range, with a spectral resolution of R=154 (Sect. <ref>). We provided an additional high-resolution XSPEC table model with an adaptive energy resolution over the 1.5 to 15 keV subrange, for fitting the microcalorimeter X-ray spectra obtained with XRISM/Resolve and Athena/X-IFU. All tables are publicly available online[<https://github.com/BertVdM/xskirtor>].

The SKIRT code can now model X-ray polarisation in AGN circumnuclear media and predict the spectro-polarimetric X-ray signal of complex 3D models, with all features of the established SKIRT framework available. SKIRT is highly optimised in terms of computational efficiency, allowing for complex 3D models to be explored in a short timeframe. Furthermore, the SKIRT code offers an unmatched geometrical flexibility for setting up simulations in full 3D, which has now become available to X-ray polarisation modelling. Finally, SKIRT can calculate polarisation maps at a high signal-to-noise (see Fig. <ref> and <ref>), which forms a powerful tool to study the geometrical effects that are encoded in spectro-polarimetric observations. SKIRT can calculate fluxes, images, spectra and polarisation maps from mm to X-ray wavelengths, and the community is warmly invited to use the code in any way they see fit.


      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. G. M. acknowledges financial support from Italian MUR under grant PNRR-M4C2-I1.1-PRIN 2022-PE9-An X-ray view of compact objects in polarized light -F53D23001230006-Finanziato dall'U.E.-NextGenerationEU. We wish to thank K. A. Arnaud for support with the XSPEC package. 

aa


