© The Authors 2025

1. Introduction

Ultraluminous X-ray sources (ULXs) were discovered in external galaxies in the 1980s with the Einstein X-ray Observatory Fabbiano (1989), Long et al. (1981). One of the most distinctive features of ULXs is the high X-ray luminosity, which can reach values of L ∼ 10⁴¹ erg s⁻¹, far exceeding the Eddington limit for a solar mass object (see e.g. Kaaret et al. 2001, for a review). From an observational standpoint, ULXs are characterized by a higher occurrence in star-forming galaxies (Swartz et al. 2004), a turnover in the spectra at energies in the 2–10 keV range, a second softer component below 1 keV, and spectral variability between different epochs (e.g. Gladstone et al. 2009; Pintore & Zampieri 2012; Sutton et al. 2013a,b; Middleton et al. 2015a,b).

The large amount of energy emitted and the X-ray variability led to believe that ULXs were binary systems in which a massive donor transfers material onto a black hole (Kaaret et al. 2017). For some time, ULXs were thought to be powered by black holes of intermediate mass (IMBHs), that is, larger than 100 M_⊙ (Colbert & Mushotzky 1999; Makishima et al. 2000; Miller et al. 2003). In fact, the Eddington luminosity for such objects exceeds the observed one, so accretion can proceed at a sub-Eddington rate. Later on, super-Eddington accretion onto stellar-mass black holes (M_BH ∼ 10–20 M_⊙, Stobbart et al. 2006; Feng & Soria 2011 and also black holes with M_BH ∼ 30–80 M_⊙, Zampieri & Roberts 2009) were also considered. If this was the case, however, emission needs to be beamed to match the observed luminosity (Fabrika & Mescheryakov 2001; King et al. 2001; Poutanen et al. 2007).

This picture was soon revolutionized in 2014, when pulsations were discovered in M82 X-2, a ULX in the M82 galaxy (Bachetti et al. 2014), supporting the view that at least some of these sources could be powered by a neutron star (NS). The discovery of the first pulsating ULX was soon followed by several others (Fürst et al. 2016; Israel et al. 2017a,b; Rodríguez Castillo et al. 2020; Carpano et al. 2018; Sathyaprakash et al. 2019; Quintin et al. 2021), highlighting the existence of an entire population of pulsating ULXs, dubbed PULXs.

Even though the association of PULXs with neutron stars appears well established, explaining how neutron stars can emit such a huge amount of energy remains an open problem, given that the Eddington luminosity for a neutron star is ≈10³⁸ erg s⁻¹. Accretion onto (highly) magnetized NSs has been the focus of many investigations, starting from the pioneering work of Basko & Sunyaev (1976). They presented a model for accretion onto magnetized neutron stars that included the formation of an accretion column on the magnetic poles of the star for the first time. This model for the column was somewhat simplistic and did not take into account the consequences of the propeller effect, induced by the huge magnetic field strength (∼10¹⁵ G), on the column. Subsequent investigations addressed the problem with an increasing level of sophistication, including the effects of a super-strong magnetic field (see Lyubarskii & Syunyaev 1988; Mushtukov et al. 2015, 2017; Brice et al. 2021). In particular, it was realized that in a strongly magnetized accretion column, the opacity is drastically suppressed for photons polarized perpendicularly to the magnetic field direction. This results in a much higher value of the Eddington luminosity, allowing for a larger flux to escape from the sides of the column. Such a scenario could indeed explain the large luminosity observed from these sources.

We present here a simplified model (torusdisk hereafter) that reproduces the X-ray thermal emission of pulsating ULXs, in terms of emission from an accretion disk and an accretion envelope, which is modeled as a torus confined by the magnetic field lines that extend up to the magnetospheric radius. To this aim, we adapted the numerical code by Taverna & Turolla (2017), which calculates the flux of photons coming from the part in view of the system as a function of both the photon energy and the star rotational phase, so to compute spectra and light curves of PULXs in different energy bands. In addition to the flux, the code also computes the polarization degree (PD) and polarization angle (PA) as functions of energy and phase.

Finally, we compare the simulations with the observational data of two PULXs, to test if the model can reproduce the properties exhibited by real sources. Firstly we focus on M51 ULX-7 (aka NGC 5194 X-7, Roberts & Warwick 2000, CXOM51 J133001.0+47134, Terashima & Wilson 2004, NGC 5194/5 ULX-7, Liu & Mirabel 2005), located in the outskirts of a young open cluster in a spiral arm of its host galaxy. It belongs to a high-mass X-ray binary (HMXB) (Abolmasov et al. 2007) with a (likely) O-B companion with mass ≳8 M_⊙ (Rodríguez Castillo et al. 2020; Abolmasov et al. 2007) and exhibits an X-ray luminosity of 6 × 10³⁹ erg s⁻¹. It was suggested that the source is powered by a neutron star accreting at a super-Eddington rate (Brightman et al. 2022). NGC 7793 P13 (P13 hereafter) is a pulsating ULX located on the southern edge of the galaxy NGC 7793 (in the Sculptor group; Fabbiano et al. 1992), the X-ray emission of which is completely dominated by this source. The ESO-VLT observations in 2008–2009 revealed that NGC 7793 P13 has a B9Ia supergiant companion with a mass between 10 and 20 M_⊙ (Motch et al. 2011; Fürst et al. 2021). From XMM-Newton and NuSTAR observations in 2016, coherent pulsations were detected (P ≈ 0.42 s); this led to identify P13 as the third known ULX powered by a neutron star (Israel et al. 2017b; Fürst et al. 2016).

The main goal of this paper is to derive information on the viewing geometry of these sources confronting the simulated light curves with the observed ones, to constrain the temperature at the inner radius of the disk, and the magnetic field strength. To reproduce the spectrum in the 0.1 − 10 keV band, three components are required: the spectrum from our model (torus and disk), a cut-off power law component to account for the non-thermal emission in the hard band (believed to come from the accretion column), and a soft blackbody component, possibly due to radiative winds produced in the intermediate/outer part of the accretion disk (Poutanen et al. 2007; Pinto et al. 2016, 2020; Middleton et al. 2022). In case spectral fitting is degenerate for different geometries, we use polarization to discriminate among them. This further shows that polarimetric measurements are indeed a powerful tool for the study of highly magnetized neutron stars in PULXs, as testified by the results obtained by the recent NASA-ASI mission IXPE (Imaging X-ray Polarimetry Explorer, Weisskopf et al. 2022).

This work is organized as follows. In Section 2 we present the model and the assumptions on which it relies. In Section 3 we provide a brief description of how the data were analyzed while Section 4 is devoted to the spectral and polarization results, which are then discussed in Section 5. Finally, our conclusions are summarized in Section 6.

2. Physical model

In this section, we introduce the physical model adopted to simulate and reproduce the X-ray spectrum and polarization properties of the radiation emitted by pulsating ULXs.

2.1. Modeling the source

The model computes the emission from an accreting magnetized neutron star (NS) in a binary system. The accreting material originates from the donor star and proceeds through a geometrically-thin accretion disk. At the Alfvén radius (r_A), where the magnetic pressure equals the ram pressure of the gas particles, the matter is funneled along the closed field lines of the magnetic field, falling towards the polar regions. This envelopes the star in a sort of cocoon that, for the sake of simplicity, we call “torus” hereafter (Mushtukov et al. 2015, 2017; Brice et al. 2021). At the high accretion rates expected to occur in PULXs, this torus is optically thick (see Fig. 1). In the following, we assume that the neutron star has a mass of 1.4 M_⊙, a radius of 10 km, and that the magnetic field has a dipolar topology with a strength (at the pole) in the range 10¹²–10¹³ G. The equation of the magnetic field lines is, in polar coordinates,

$$\begin{array}{c}\hfill r=r_{\mathrm{A}}{sin}^2\theta ,\end{array}$$(1)

Fig. 1. A representation of the system. A 3D view on the top and a meridional cut on the bottom. The accretion disk and torus are shown in orange and blue respectively, the spin axis Ω and the magnetic dipole axis bdip are also indicated (in gray the radiation emitted from the accretion column).

where θ is the magnetic colatitude, and r_A, for the sake of simplicity, here coincides with the maximum radius of the magnetosphere, which is computed with the equation (1) in Brice et al. (2023):

$$\begin{array}{c}\hfill R_{\mathrm{m}}=7\times {10}^7\mathrm{\Lambda }{M}^{1/7}{R}_6^{10/7}{B}_{\mathrm{d},12}^{4/7}{L}_{39}^{-2/7}\mathrm{cm},\end{array}$$(2)

where Λ is a dimensionless parameter that depends on the mode of accretion (a typical value of 0.5 for accretion via thin disc), B_d, 12 is the surface dipole field strength at the magnetic poles (expressed in units of 10¹² G), and L₃₉ is the accretion luminosity (in units of 10³⁹ erg s⁻¹). The equation (1) is just an approximation since the exact topology of the magnetic field near the surface of the NS is unknown. Pulsating ULXs, in particular, may have a multipolar magnetic field structure, as proposed by Israel et al. (2017a) to accommodate for the large field required to rise the Eddington limit enough and, at the same time, prevent the star from entering the propeller stage. However, higher-order multipoles quickly decay moving away from the surface, so our assumption is quite justified if r_A is tens stellar radii and makes the computations much simpler.

Both the disk and the torus are assumed to be in Local Thermal Equilibrium (LTE). The disk temperature follows the thin-disc profile (see Shakura & Sunyaev 1973)

$$\begin{array}{c}\hfill T(r)=T_{\mathrm{in}}{\left(\frac{R_{\mathrm{in}}}{r}\right)}^{3/4},\end{array}$$(3)

where T_in is the temperature at the inner radius of the disk, R_in, that we take to coincide with the radius at which the torus and the disk intersect. The temperature profile of the torus is that presented in Brice et al. (2023), it varies with the magnetic colatitude and follows the expression:

$$\begin{array}{c}\hfill T_{\mathrm{out},\mathrm{torus}}=T_{\mathrm{in},\mathrm{torus}}{\tau }^{-1/4},\end{array}$$(4)

where T_out, torus is the temperature at the torus surface, T_in, torus the temperature at the torus inner boundary, and τ > 1 the torus optical depth (see Figure 1 in Brice et al. 2023). As reported in Brice et al. (2023), the local optical depth depends on the dynamics of the accreting particles: the non-uniform velocity field inside the envelope introduces a meridional variation of T_out. Because of our assumption of LTE, each annulus of both the disc and the torus emits blackbody radiation at the local temperature; no effects due to the reprocessing of thermal radiation by material on top of the disc and/or torus are considered. In our simplified modeling, we do not account for radiation coming from the sides of the accretion columns which form above the magnetic poles if the mass accretion rate is high enough. Emission from the columns gives rise to a non-thermal tail in the X-ray spectrum (Walton 2015; Walton et al. 2018a; Brightman et al. 2022). In order to reproduce the observed spectrum over the entire range, we added to the model a phenomenological cut-off power law, often used to describe this type of non-thermal emission (see e.g. Walton et al. 2018b).

2.2. Polarization

Radiation emitted by strongly magnetized neutron stars is expected to be highly polarized in two normal modes, the ordinary (O) and extraordinary (X) ones (e.g. Gnedin & Pavlov 1974; Ho & Lai 2003; Lai et al. 2010). The former is characterized by an electric field oscillating in the plane of the propagation vector k and the local magnetic field B, while in the latter the electric field oscillates perpendicularly to both these vectors. The cross-sections of X-mode photons are greatly reduced below the electron cyclotron energy, E_c, e ≈ 11.6(B/10¹² G) keV, by a factor of ∼(B/B_Q)² with respect to the unmagnetized case, where B_Q ≃ 4.414 × 10¹³ G is the critical field (Mészáros 1992; Harding & Lai 2006). This makes the medium optically thin for X photons (Herold 1979; Ventura 1979) and leads to the release of more radiation inside the accretion column, explaining the super-Eddington luminosity in the X-rays.

Moreover, strong magnetic fields are also expected to modify the optical properties of the medium in which radiation propagates. In particular, vacuum birefringence (Heisenberg & Euler 1936) affects photons propagating in a strongly magnetized vacuum, forcing the polarization modes to remain unchanged within a region close to the star surface where the magnetic field is strong enough (Heyl & Shaviv 2002; Heyl et al. 2003; Fernández & Davis 2011; Taverna et al. 2014). Inside the so-called adiabatic region, the scale length ℓ_A along which the photon electric field varies turns out to be much smaller than the scale length ℓ_B that characterizes the variation of the star magnetic field along the photon trajectory. As a result, the photon electric field adapts instantaneously to the local magnetic field direction. On the other hand, far from the star, it is ℓ_A ≫ ℓ_B, so that the electric field is frozen. The boundary between these two regimes is given by ℓ_A ≃ ℓ_B, which is called adiabatic (or polarization-limiting) radius (Heyl et al. 2003; Taverna et al. 2015):

$$\begin{array}{c}\hfill \frac{r_{\mathrm{pl}}}{R_{\mathrm{NS}}}\sim 5{\left(\frac{E}{1\mathrm{keV}}\right)}^{1/5}{\left(\frac{R_{\mathrm{NS}}}{10\mathrm{km}}\right)}^{1/5}{\left(\frac{B_{\mathrm{p}}}{{10}^{11}\mathrm{G}}\right)}^{2/5},\end{array}$$(5)

where E is the photon energy, B_p the polar magnetic field strength, and R_NS the NS radius. The actual polarization state should be determined by solving the wave equation. However, here we adopt the approximation already discussed in Taverna et al. (2015) and assume adiabatic evolution for r < r_pl while the electric field is frozen for r > r_pl. For the magnetic field strength (∼10¹³ G) and the photon energies considered here (0.1–1 keV), the Alfvén radius (r_A) is larger than the adiabatic radius (r_A ∼ 50R_NS, while r_pl ∼ 36R_NS), so we expect that a large part of the torus lies outside the adiabatic region. Hence, we assume radiation coming from this zone of the torus to be unpolarized¹. Radiation coming from inside r_pl carries a non-zero intrinsic polarization degree,

$$\begin{array}{c}\hfill {\mathrm{\Pi }}_{\mathrm{L}}=\left|\frac{N_{\mathrm{X}}-N_{\mathrm{O}}}{N_{\mathrm{X}}+N_{\mathrm{O}}}\right|,\end{array}$$(6)

where N_X (N_O) is the fraction of X (O) photons. Following the same argument, we assume that the polarization of radiation coming from the disk is negligible since the disk starts outside r_pl where B is quite small².

It is convenient to express the polarization observables, the polarization degree (PD) and polarization angle (PA), in terms of the Stokes parameters ℐ, 𝒬, and 𝒰 (Rybicki & Lightman 1991). In a reference frame with the z axis along the unit wavevector k, and the y axis in the (k, B) plane, the Stokes parameters for X and O photons are

$$\begin{array}{c}\hfill {\left(\begin{array}{c}\mathcal{I}\\ \mathcal{Q}\\ \mathcal{U}\end{array}\right)}_X=\left(\begin{array}{c}1\\ 1\\ 0\end{array}\right){\left(\begin{array}{c}\mathcal{I}\\ \mathcal{Q}\\ \mathcal{U}\end{array}\right)}_O=\left(\begin{array}{c}1\\ -1\\ 0\end{array}\right).\end{array}$$(7)

The last Stokes parameter, 𝒱, which describes the circular polarization, is not considered here. The Stokes parameters depend on the reference frame of each photon (x_i, y_i, z_i), which is fixed by the local direction of the (dipole) magnetic field that changes across the emitting surface. To sum the Stokes parameter we must refer them to the same fixed reference frame, that we choose to be that of the polarimeter. Each photon reference frame is rotated (around z_i) with respect to the polarimeter frame by an angle α_i. Eventually, for a single photon the Stokes parameters I_i, Q_i, U_i are defined as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill I_i& ={\mathcal{I}}_i\hfill \\ \hfill Q_i& ={\mathcal{Q}}_{i}cos(2{\alpha }_i)+{\mathcal{U}}_{i}sin(2{\alpha }_i)\hfill \\ \hfill U_i& ={\mathcal{U}}_{i}cos(2{\alpha }_i)-{\mathcal{Q}}_{i}sin(2{\alpha }_i).\hfill \end{array}\end{array}$$(8)

The observed PD and PA are given by Rybicki & Lightman (1991)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathrm{PD}& =\frac{\sqrt{Q^2+U^2}}{I}\hfill \\ \hfill \mathrm{PA}& =\frac{1}{2}arctan\left(\frac{U}{Q}\right),\hfill \end{array}\end{array}$$(9)

where the Stokes parameters Q and U refer to the “total” radiation, that is, the sum of the Q_i and U_i of the collected photons. Substituting Equations (7) and (8) inside Equations (9) one obtains

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \begin{array}{cc}\hfill \mathrm{PD}& =\frac{1}{N}[N+2{\displaystyle \sum _i}{\displaystyle \sum _{k>1}}cos(2{\alpha }_i-2{\alpha }_k)\hfill \\ \hfill & +2{\displaystyle \sum _j}{\displaystyle \sum _{r>j}}cos(2{\alpha }_j-2{\alpha }_r)\hfill \\ \hfill & -2{\displaystyle \sum _i}{\displaystyle \sum _j}cos(2{\alpha }_i-2{\alpha }_j){]}^{1/2}\hfill \\ \hfill \mathrm{PA}& =\frac{1}{2}arctan-\frac{{\sum }_{i}sin(2{\alpha }_i)-{\displaystyle {\sum }_j}sin(2{\alpha }_j)}{{\sum }_{i}cos(2{\alpha }_i)-{\displaystyle {\sum }_j}cos(2{\alpha }_j)},\hfill \end{array}\end{array}\end{array}$$(10)

where i, k = 1, …, N_X and r, j = 1, …, N_O. From Equation (10) one can see that the polarization degree measured at infinity is in general lower than that at the emission (Equation 6), because of the presence of the geometrical factors sin(2α) and cos(2α). If the magnetic field topology is particularly tangled (as it occurs close to the star surface), the α angles of the different photons span almost over the entire 0–2π range and the depolarization effect is larger. On the other hand, for a more uniform magnetic field (like that at r_pl, far from the surface), the different α angles attain much more similar values and radiation is much less depolarized (Taverna et al. 2015). Hence, when vacuum birefringence is accounted for and the Stokes parameters are computed at r_pl, the polarization pattern at the surface is likely preserved also at the observer.

2.3. Numerical implementation

To compute the spectral and polarization properties of the emitted radiation, we used the ray-tracer code discussed in Taverna et al. (2015), but adapted to our model, in particular to a different emitting surface (torus and disk).

2.3.1. Visibility of the source

The code starts by computing the visible part of the source assuming a specific viewing geometry. We consider the surface, comprising the disk and the torus, as a collection of small emitting patches. The grid is made by 50 bins in θ, ϕ, and the radial distance r, while for both the energy and rotational phase we take 30 bins. The phase is measured from the projection of the line of sight (LOS) onto the plane perpendicular to the spin axis.

In the LOS reference frame (with the x-axis in the plane of l and Ω, see Fig. 2), the code selects the points that are in view for the chosen angles χ and ξ, which represent the orientations of the neutron star spin axis and the magnetic axis with respect to the LOS, respectively.

Fig. 2. LOS reference frame. The Z-axis is aligned to l, the LOS unit vector, η is the angle between the star magnetic axis bdip and the LOS, while ζ is the corresponding azimuth. The LOS and magnetic axis inclinations χ and ξ with respect to the spin axis Ω are also shown.

The code also takes into account the orientation of the disk relative to the NS spin axis. The following steps are performed: selecting the visible part of the torus, accounting for its self-shadowing as well as the shadow cast by the disk, and selecting the visible part of the disk not shadowed by the torus.

To achieve this, the code traces a straight line parallel to the LOS and identifies the intersections between this line and the emitting surface (torus and/or disk). To understand if each intersection is in view, we evaluate the quantity

$$\begin{array}{c}\hfill \mathit{n}·\mathit{l}>0,\end{array}$$(11)

where n represents the surface normal unit vector at the intersection, and l is the unit vector of the LOS. If multiple intersections satisfy Equation (11), we calculate the z-coordinate of each point and consider in view the one with the largest value of the z-coordinate.

Figure 3 shows the emitting regions of the source for a given viewing geometry, in the LOS reference frame. Although a full treatment requires to take into account general relativistic (GR) effects, such as gravitational redshift and ray-bending, we expect them to be negligible since the emitting regions (torus and disc) lie in general far from the star. For these reasons, we do not consider GR corrections in the present simulations.

Fig. 3. The selection of the part in view for the torus, in panels (a) and (b), while considering also the disk, in panels (c) and (d). In panels (a) and (c) the source is shown exactly along the LOS direction, while in panels (b) and (d) the reference frame is rotated, for ease of visualization. The points in view are marked in cyan for the torus and in orange for the disk, while those not in view are marked in yellow for the torus and white for the disk. In this case, we used χ = 45°, ξ = 1°, Rmax = 70 RNS, and Rdisk = 200 RNS.

2.3.2. The flux

Once a point within the emission region of the torus is tagged as in view, we calculate its temperature. As reported above, the temperature of the torus varies with the magnetic colatitude. The temperature of each emission point is obtained by interpolating over the grid discussed in Brice et al. (2023). Instead, for a visible point on the disk we derive the temperature from the Shakura-Sunyaev profile (Equation 3), once the inner radius R_in and the inner temperature T_in of the disk are fixed (Section 4.2). The locally emitted flux is calculated from the blackbody distribution at the local temperature. We then sum all the contributions coming from different points in view using the same approach as discussed in Taverna & Turolla (2017, see their Section 3.2), obtaining as output a matrix with the total observed flux as a function of energy and phase.

2.3.3. Polarization observables

For each point in view inside r_pl we set N_X (N_O = 1 − N_X) and compute the rotated Stokes parameters of the photons at the adiabatic radius (Equation (5) in Taverna et al. 2015). As described above for the total flux, we sum the Stokes parameter fluxes over the part in view of the emitting region (at each energy and phase), obtaining the total Q and U fluxes for both X and O photons. Hence, the code calculates PD and PA following Equation (9). For the points in view that are at a distance ≳r_pl from the star we do essentially the same. However, since, as described above, outside the adiabatic region the photon electric field direction is frozen with respect to the magnetic field one, we anyway expect to observe a much lower polarization degree for radiation emitted from these points. For this reason, and for the sake of simplicity, for photons coming from outside r_pl we fixed directly N_X = 0.5 (considering the 50% of radiation to be X-mode photons, and the other 50% to be O-mode photons, with no privileged polarization mode), which, from Equation (6) gives Π_L = 0. So, in this work, we analyzed the variation of PA and PD with the energy and phase varying N_X only for photons coming from the regions of the torus closer to the star surface (and so emitted from inside r_pl).

3. Data reduction

We analyzed a sample of XMM-Newton and NuSTAR observations of M51 ULX-7 and NGC 7793 P13, two PULXs taken as benchmarks for our study. We used three XMM-Newton observations for M51 ULX-7 (Obs.ID: 0824450901, 0830191501 and 0830191601), carried out between 2018 May and June, and two for NGC 7793 P13 (Obs.ID: 0693760401 and 0748390901), carried out between 2013 November and 2014 December. We extracted the EPIC-pn data using the SAS v.14.0. package³. Spectra were obtained by selecting events with the PATTERN ≤ 4 for EPIC-pn (single- and double-pixel events) and setting ‘FLAG = 0’ to ignore bad pixels and events coming from the CCD edges. Epochs of high background were also removed from the analysis. Source and background events were extracted from circular regions with radii of 30″ and 60″, respectively. The spectra were rebinned to have at least 25 counts per energy bin.

4. Results

We started by comparing the light curves of the XMM-Newton observations with the simulated ones to constrain the geometry of view, which significantly affects the flux modulation. In particular, we compared the observed and simulated pulsed fractions (PFs), defined as

$$\begin{array}{c}\hfill \mathrm{PF}=\frac{\max (\mathrm{Flux})-\min (\mathrm{Flux})}{\max (\mathrm{Flux})+\min (\mathrm{Flux})}.\end{array}$$(12)

After constraining the viewing geometry, we compared the spectra determining the temperature at the inner radius of the accretion disc (T_in), the torus temperature range, and the magnetic field strength. For the comparison of the pulse profiles, we considered six possible viewing geometries while, for the comparison of the spectra, we chose 16 values of T_in, selected on the basis of values found in the literature, for each different value of the magnetic field strength. We then analyzed the polarization degree and polarization angle for each source using the best-fit viewing geometry, T_in, and magnetic field strength.

4.1. Comparison of the light curves

We made a direct comparison of the light curves, considering the 2–3 keV and 3–4 keV energy bands since we expect that at these energies the torusdisk model (in particular the emission of the torus) dominates. We compared each observed profile with six simulations having different geometries of view, each differing only for the angle χ (ranging from 10° to 60°). The value of ξ was fixed at 10°, because for values of 20°, or higher, and for values lower than 10°, the PFs were not compatible with that observed. So, for the sake of simplicity, we fixed ξ at 10°, also to represent a general non-aligned magnetic rotator case. χ angles beyond 60° were not considered, because we noted in the simulations that the PF decreases if χ > 60°.

4.1.1. M51 ULX-7

To perform simulations we fixed some input parameters: the magnetic field strength B, the Alfvén radius r_A, the torus temperature, and the disk radius R_disk. M51 ULX-7 has a surface magnetic field strength between 8 × 10¹¹ G and 10¹³ G (values derived by Rodríguez Castillo et al. 2020, from the P-Ṗ relation), so we set a lower and upper limit, B = 10¹² G and B = 8 × 10¹² G, respectively. The Alfvén radius (r_A, which coincides with the inner radius of the disk R_in) and the temperature of the torus were calculated according to B. The outer radius of the disk was set to 200R_NS, as beyond this radius the emitted flux falls below 0.1 keV, outside the range considered here. We report all the input fixed parameters in Table 1.

For M51 ULX-7, we found a lower and upper limit for the geometry of view, χ = 20° ξ = 10° and χ = 60° ξ = 10°, respectively (see Figure 4). In Table 2 we report the values of the observed PFs with those of the best-fit viewing geometries, across the two different energy bands. The observed PFs show general agreement with the simulated ones within uncertainties.

Fig. 4. Comparison of the three XMM-Newton light curves of M51 ULX-7 with simulations. From top to bottom geometry of view with χ = 60° ξ = 10° (coral), with χ = 20° ξ = 10° (orange), and with χ = 30° ξ = 10° (green). The column on the left refers to the energy band 2–3 keV, while that on the right to the energy band 3–4 keV.

In Figure 4 we show the comparisons.

4.1.2. NGC 7793 P13

We did the same for NGC 7793 P13, considering that the observed surface magnetic dipole field strength is around 2 × 10¹² G (Fürst et al. 2016; Israel et al. 2017b, still obtained from P and Ṗ). We then set the dipole magnetic field strength to B = 4 × 10¹² G considering the value at the poles, and computed the Alfvén radius of the torus, the torus temperatures (see Table 3), and set the disk outer radius to R_disk = 200 R_NS.

This time we found much less variation for the angle χ, identifying χ = 60° ξ = 10° and χ = 40° ξ = 10° as the upper and lower limit, respectively. In Figure 5 the light curves and in Table 4 the pulsed fractions, across the two energy bands, are confronted with the observed ones.

Fig. 5. Comparison of the two XMM-Newton light curves of NGC 7793 P13 with simulations. From left to right the geometry of view with χ = 60° ξ = 10° (coral) and with χ = 40° ξ = 10° (red). The column on the right refers to the energy band 2–3 keV, while that on the left refers to the energy band 3–4 keV.

4.2. Comparison of the spectra

4.2.1. M51 ULX-7

In Section 4.1.1 we derived a lower and upper limit for the geometry of view of the source, χ = 20°, ξ = 10° and χ = 60°, ξ = 10°, respectively. We now perform spectral comparisons to constrain T_in for both the geometries and the magnetic field strengths. From the literature, the temperature of the soft spectral component (which is related to T_in) varies between 0.33 and 0.5 keV (Rodríguez Castillo et al. 2020; Brightman et al. 2022; Vasilopoulos et al. 2020). We then produced 16 spectral simulations, each with a different temperature value T_in, which ranges between 0.25 and 0.5 keV.

4.2.2. M51 ULX-7 – B = 10¹² G

We reproduced the spectrum with the superposition of two thermal components, the XSPEC BBODY model (see Arnaud 1996 and our simulated multicolor blackbody, similarly as in Rodríguez Castillo et al. 2020). We considered the interstellar absorption through the WABS model (Morrison & McCammon 1983) leaving the column density as a free parameter of the fit. We carried out the fitting procedure outside XSPEC, using a Python script where we implemented the analytical expressions of BBODY. This is done to include the TORUSDISK contribution (i.e. the output of the ray-tracer code). For the same reason, we use WABS and not TBABS, as an analytical expression is available for the former. We fixed the emission temperature of the BBODY at 0.25 keV, to ensure the convergence of the fit⁴. Furthermore, we stress that our purpose is to test the reliability of the TORUSDISK model by varying the contribution of the internal disk temperature T_in for a given viewing geometry, while accounting for the presence of additional components. In Figure 6 we show the spectral fits for both the geometries of view of the first observation (the others provide similar results). In Table 5 we report the best-fit parameters of all the observations. The best value for the temperature T_in turns out to be 0.5 keV, in agreement with Rodríguez Castillo et al. (2020) within the uncertainties. Despite the superposition of these two thermal components can fit the data, it also leaves non negligible residuals between 1.5–3 keV.

Fig. 6. X-ray spectrum of M51 ULX-7 obtained with XMM-Newton (observation: 0824450901). The spectrum is fitted with a multicomponent model. The panel on the left refers to the viewing geometry χ = 20° ξ = 10°, that on the right to χ = 60° ξ = 10°. One component is phenomenological (BBODY), and one is calculated with the model (TORUSDISK).

4.2.3. M51 ULX-7 – B = 8 × 10¹² G

In this case, we reproduced the spectrum by the superposition of a soft X-ray thermal component plus a hard X-ray non-thermal one. We fitted the former with the sum of two thermal spectral components, the XSPEC BBODY model and the TORUSDISK multicolor blackbody, while the latter with a cut-off power law (CUTOFFPL in XSPEC), accounting for the interstellar absorption through the WABS model. We carried out the fitting procedure, as before, outside XSPEC, implementing, in this case, also the analytical expressions of CUTOFFPL. Still to ensure the convergence of the fit, we fixed the temperature of the BBODY and also the column density, N_H = 9 × 10¹⁹  cm⁻²⁵.

We used the superposition of a blackbody, a multicolor blackbody and a power law component as done by Brightman et al. (2022). The low-energy blackbody component, in fact, is used to describe the emission due to radiative winds that can originate from the accretion disk at super-Eddington accretion rates, (see e.g. Poutanen et al. 2007; Pinto et al. 2016, 2020). In Figure 7 we show the spectra for the three XMM-Newton observations. In Table 6 we report the best-fit values of the parameters. The value of T_in turns out to be ∼0.3 keV, compatible with those on the literature, for both the viewing geometries, which also show indistinguishable spectral results. Considering the results with B = 8 × 10¹² G, in the following, we took into account only the case with this value of the magnetic field strength.

Fig. 7. X-ray spectra of M51 ULX-7 obtained with XMM-Newton (observations: 0824450901, 0830191501, 0830191601). The spectrum is fitted with a multicomponent model. The panels on the left refer to the viewing geometry χ = 20° ξ = 10°, while those on the right to χ = 60° ξ = 10°. Two components are phenomenological (BBODY, CUTOFFPL), and one is computed with the model (TORUSDISK).

4.2.4. NGC 7793 P13

We adopted the same approach for the spectral fits on NGC 7793 P13, using the geometries χ = 60° ξ = 10° and χ = 40° ξ = 10° derived form the pulse profiles analysis. From the X-ray spectral analyses of Walton et al. (2018a), the spectrum of this source can be fitted with a model made by three components, two of which are thermal and dominate below 10 keV. They found that the temperature of the coolest thermal component spans from ∼0.3 to 0.5 keV, and that of the hottest thermal component from ∼1 to 1.5 keV. We used the same multicomponent model described above (bbody + torusdisk model + cutoff power law), correcting for absorption (WABS) with a fixed column density, N_H = ∼3 × 10²⁰  cm⁻², and temperature of the BBODY component, 0.25 keV. We then varied T_in between 0.25–0.5 keV in 16 bins. We show the results in Figure 8 and we report the values of the parameters in Table 7. We find T_in ∼ 0.3 keV, in agreement with the literature, and also in this case the spectra of the two geometries of view do not show differences. The results also suggest that B = 4 × 10¹² G can reproduce the spectrum well.

Fig. 8. X-ray spectrum of NGC 7793 P13 obtained with XMM-Newton (observations: 0693760401 on the top, 0748390901 on the bottom). The spectrum is fitted with a multicomponent model. Two components are phenomenological (BBODY, CUTOFFPL), and one is computed with the model (TORUSDISK). Panels on the top (bottom) from left to right are comparisons with geometries χ = 40° and χ = 60°, respectively.

4.3. Polarization

In simulating the polarization properties of the two PULXs, we considered fields with strengths of 4 × 10¹² G (NGC 7793 P13) and 8 × 10¹² G (M51 ULX-7) at the poles. For photons escaping from the adiabatic region, we fixed the intrinsic polarization fraction, Π_L, at 60% in the X-mode. We did not computed this value self-consistently since this would have required to solve the full radiation transport problem, which is definitely outside the scope of the present work. Nevertheless, in the presence of strong magnetic fields (≳10¹³ G) X-mode photons are expected to dominate the emission, so our assumption is not unreasonable. The effects of changing Π_L will be discussed further on. We then analyzed the variation of the polarization degree and polarization angle with the geometry of view. In the case of M51 ULX-7 we used χ = 60°, 40° and 20°, while for NGC 7793 P13 χ = 60° and 40°.

4.3.1. M51 ULX-7

Fig. 9. PD and PA variation with energy and phase for different geometries of view. From left to right χ = 60°, 40°, and 20°.

We show in Figure 9 the energy- and phase-dependent PD, for three different viewing geometries. In each panel, the polarization degree increases with the energy and shows an oscillating behavior with the phase, that changes going from χ = 20° to χ = 60°. The former is a consequence of the strong magnetic field considered in the model, while the latter is a consequence of the rotation of the source combined with the geometry of view.

Fixing energy and phase but varying the angle χ from 20° to 60° the value of the polarization degree increases. This effect could be caused by the geometrical depolarization. The polarization angle oscillates with the phase remaining almost constant with the energy. As for the polarization degree, the variation with the phase of the polarization angle depends on the viewing geometry.

4.3.2. NGC 7793 P13

In Figure 10 we show the polarization results for NGC 7793 P13.

Fig. 10. PD and PA variation with energy and phase for different geometries of view. From left to right χ = 60° and 40°.

The maximum polarization degree is lower than that for M51 ULX-7. Fixing energy and phase, it is higher for the geometry of view with χ = 40° than that with χ = 60°. As for M51 ULX-7, the polarization degree oscillates with the phase and increases with the energy, while the polarization angle oscillates with the phase remaining quite constant with the energy for low values of χ.

4.3.3. Polarization results with different Π_L

We previously mentioned that the value of Π_L is not computed, but it is fixed considering the magnetic field strength of the source. We report here the results for the polarization degree varying Π_L and fixing the geometry of view with χ = 60° and ξ = 10°, for both M51 ULX-7 and NGC 7793 P13. The results show that a variation in the intrinsic polarization degree only leads to an overall variation in the value of the PD, without changing its behavior with energy and phase. This reinforces the hypothesis that, even if the intrinsic polarization value assumed in the model may be different from the one computed self-consistently, we can still distinguish the geometry of view through the energy and phase variations of PD, since they do not depend on the value of Π_L, but only on the geometry of view and magnetic field of the source.

Fig. 11. Polarization degree variation for different values of ΠL, for M51 ULX-7 with B = 8 × 1012 G and geometry of view χ = 60° ξ = 10°.

Fig. 12. Polarization degree variation for different values of ΠL, for NGC 7793 P13 with B = 4 × 1012 G and geometry of view χ = 60° ξ = 10°.

5. Discussion

This paper aims to constrain the geometry of view, the magnetic field strength, and thermal properties of pulsating ULXs, reproducing the thermal radiation emitted by a highly magnetized neutron star accreting at a super-Eddington rate and taking into account the emission from an accretion disk and an optically thick envelope surrounding the magnetosphere. The model reproduces the light curves and the spectra of PULXs, for the latter by adding one or two phenomenological components that are not self-consistently included in the model: a blackbody component to model the excess of flux at low energies, probably produced by optically thick winds from an intermediate region of the disk; a cut-off power law produced by an accretion column. In addition to the spectral properties, we also incorporated the polarization observables of the emitted radiation, focusing our analyses on the polarization degree. We discuss our results in the following.

5.1. M51 ULX-7

5.1.1. Light curves and spectra

We started comparing the light curves because the variation of the flux with the phase, and the corresponding pulsed fraction, are highly dependent on the viewing geometry, allowing us to constrain it. We considered light curves in two energy bands, 2–3 keV and 3–4 keV, where we expect our model, in particular the torus emission, to dominate. Being the emission from the accretion column, which contributes to the PF at high energies, not included in our model, it would be impossible to compute the PF self-consistently. For M51 ULX-7 we did not find a unique solution for the viewing geometry but we were able to derive a lower and upper limit, χ = 20°, ξ = 10° and χ = 60°, ξ = 10°, respectively. This is not surprising since the source is very variable and its PF does not increase monotonically with the energy, as shown in Figure 3 in Brightman et al. (2022). The reasons behind the source variability are not well known. They could be related to the fact that the source entered a propeller stage, or because of a variation in the mass transfer rate, which is unstable.

We used both the geometries to perform the spectral comparison, adding to the torusdiks multicolor blackbody the phenomenological components to fit the entire spectrum. We considered two representative values for the magnetic field, B = 10¹² G and 8 × 10¹² G, respectively, which brackets the true value. For each of them, we computed the Alfvén radius r_A and the torus temperature. We also fixed the disk outer radius, R_disk = 200 R_NS, and varied the temperature at the inner radius of the disk between 0.25 and 0.5 keV.

For a magnetic field B = 10¹² G we were able to fit the spectrum with two thermal components: a blackbody at lower energies and our torusdik model at higher energies (Figure 6), similar to what was done by Rodríguez Castillo et al. (2020). In the fitting procedure, we left free to vary the normalization of the torusdisk model and the column density N_H, while we fixed at 0.25 keV the BBODY temperature, and we let the temperature T_in to vary between 0.25 and 0.5 keV, but fixing it for each simulated spectrum. We did not fix the parameters of the phenomenological BBODY component to the values inferred by Rodríguez Castillo et al. (2020), because our multicolor blackbody is different. The two thermal component models of Rodríguez Castillo et al. (2020) fit well the spectrum, with a higher column density of 5.9 × 10²⁰ cm⁻², and a higher value for the temperature and radius of the hotter thermal component (1.5 keV and 97 R_NS, respectively), while T_in is compatible within the uncertainties with our value. On the other hand, our physically motivated model (blackbody + torusdisk) does not provide a satisfactory fit to the data (Figure 6) because of residuals between ∼1.5 and 3.0 keV. The reason is that in our torusdisk multicolor blackbody, the energy range of the torus temperature varies depending on the choice of the magnetic field strength and cannot be adjusted in the fit. In addition, Rodríguez Castillo et al. (2020) also adopted the absorption model TBABS, an updated version of WABS, and this can explain the higher N_H.

For a magnetic field B = 8 × 10¹² G, we fitted the data with the superposition of three components: a phenomenological blackbody, the torusdisk model, and a phenomenological cutoff power law, which represents the non-thermal emission from the accretion column. With the three component model, the spectrum is reproduced with great accuracy (Figure 7). A similar 3-component spectral analysis was also performed by Brightman et al. (2022). They fitted the data using the same blackbody and cutoff power law components, and a phenomenological multicolor blackbody at intermediate energies, well reproducing the whole spectrum. For the same reasons repeated above we did not use their best-fit parameters for the blackbody and the cutoff power law, but we left them free to vary (apart from those that we fixed to ensure the convergence of the fit: the BBODY temperature and the total column density). The T_in is ∼0.36 keV, compatible with our value, while the column density N_H is lower than that found by Rodríguez Castillo et al. (2020), N_H ∼ 3.3 × 10²⁰ ± 0.7 cm⁻², and closer to our value (the residual difference being again possibly related to the different absorption model). For the cutoff power law parameters, Brightman et al. (2022) fixed them using the values found by Walton et al. (2018a), E_cut(here β)∼8.1 keV and Γ(here α)∼0.8, that are different from our results.

From spectral comparison we obtained good constraints on the temperature T_in and on the magnetic field strength B, but no additional constraints on the geometry of view. The latter can be further investigated through polarization analyses.

5.1.2. Polarization

As shown in Sect. 4.3.1 the polarization degree increases with the energy (Figure 9). This is a consequence of the strong magnetic field considered. For B = 8 × 10¹² G, high energy photons are emitted from the torus inside the adiabatic region. This explains the oscillating behavior of the polarization degree with the phase, which is simply a consequence of the rotation of the source combined with the geometry of view, that changes the emitting regions in view from which the polarized radiation comes. For regions closer to the equator PD is low because photons have lower energies (the torus temperature is lower) and are emitted outside the adiabatic region.

Fixing the energy and the phase, the polarization degree increases with the angle χ (Fig. 9). This is because of the geometrical depolarization effect. To understand this point, let us consider two limiting cases: a pole-on configuration (χ = ξ) and an equator-on configuration (χ − ξ = 90°), with respect to the magnetic axis. In the first case, there is no preferential direction of the magnetic field vector, since the field assumes a radial configuration. Therefore photons escaping from those regions (linearly polarized parallel to or perpendicular to the plane of the vectors k and B) will have an electric field with a radial configuration, which means to have opposite polarization contribution, so a null total polarization. In the second case, there is a preferential direction of the magnetic field vector, so summing the polarization contributions we have a non-zero polarization degree. For the χ = 20°, ξ = 10° geometry, we are closer to a pole-on configuration, while increasing χ up to 60° we move towards an equator-on configuration.

It is also worth noting that even if we cannot distinguish the viewing geometry from the X-ray spectrum, polarization measurements may allow us to do it. Indeed, the behavior of the polarization degree remains the same even changing the intrinsic polarization fraction Π_L.

The polarization angle shows an oscillating behavior with the phase, going from 90° to 0° or 180°, as expected from a rotating vector model, while remaining almost constant with the energy. The oscillations are caused by the geometry of view combined with the rotation of the star. Because the part in view of the emitting components varies from regions emitting polarized radiation (X-mode photons), where we expect the polarization angle to oscillate around 90°, to regions emitting mainly unpolarized radiation, that give a PA of 0° or 180°, if we recall how PA is defined (Equation (9)), and that the values of U and Q, in those regions, are close to zero.

5.2. NGC 7793 P13

5.2.1. Light curves and spectra

We repeated the previous analysis on NGC 7793 P13, starting from the light curve comparisons. We found two limiting viewing geometries (this time with a smaller range of χ): χ = 40°, ξ = 10°, and χ = 60°, ξ = 10°. The pulsed fraction of NGC 7793 P13 presents a more regular behavior, increasing quasi-monotonically with energy (Israel et al. 2017b). We used both geometries to perform the spectral fits, considering a magnetic field strength of 4 × 10¹² G (computing the torus temperature, the Alfvén radius, Table 3, and range of T_in as before), and setting the disk outer radius R_disk = 200 R_NS. The spectrum is well reproduced by the superposition of a blackbody, the torusdisk model, and a cutoff power law. Israel et al. (2017b) performed a phenomenological fit using an absorbed power law with a high energy cut-off and a multicolor blackbody (PHABS[HIGHECUT*(POWERLAW+BBODYRAD)]), obtaining the following parameters: N_H = 9.6 × 10²⁰ cm⁻², Γ(here α)∼1.2, E_cut(here β)∼6 keV, kT_BB ∼ 0.2 keV. Despite some differences in the fitting parameters, both models agree on the type of emission at high energies, well described by a power law, probably produced by an accretion column, while the emission in the soft X-ray energy band is thermal and well described by a multicolor blackbody.

5.2.2. Polarization

As for M51 ULX-7, also in NGC 7793 P13 for a given geometry of view the polarization degree increases with energy, because of the strong magnetic field, and it oscillates with phase. The maximum polarization degree is around 4.7%, lower than that of M51 ULX-7 (∼19%), using the lower magnetic field strength, B = 4 × 10¹² G. For such magnetic field we expect the radiation to be poorly polarized (in the X-mode), expecting less variation of PD between more energetic and less energetic radiation (Figure 10, top right). We observe oscillations with the phase still related to the rotation of the source combined with the geometry of view. The reason why this time the polarization degree decreases with χ depends on the different magnetic field strength used, on the geometry of view, and on how the geometrical depolarization effect dominates over the vacuum birefringence and vice versa. Indeed, despite the favorable viewing angle (χ = 60°), because of the depolarization effect mentioned above, the lower magnetic field strength reduces the polarized radiation. As a result, in regions near the equator the radiation is so weakly polarized that, even summing all the polarization contributions, we cannot achieve a higher polarization degree (Figure 10, top left).

Also for NGC 7793 P13, the behavior of PD may allow us to distinguish the geometry of view if observed polarization measurements are available.

The polarization angle oscillates with phase and remains constant with energy, but switching drastically between the highest (∼180°) and lowest value (∼1°) at lower energies for the geometry of view χ = 60° ξ = 10°. The reason is probably that the radiation is mainly unpolarized.

6. Conclusions

In this work, we presented a simplified model that reproduces the X-ray thermal emission of pulsating ultraluminous X-ray sources, in terms of emission from an accretion disk and an accretion envelope modeled as a torus confined by the magnetic field lines and extending up to the magnetospheric radius. We used and modified the ray-tracer code discussed in Taverna et al. (2015), to reproduce spectra and light curves of PULXs in two different energy bands. The model also simulates the polarization degree and angle of the radiation emitted as a function of energy and phase. We tested the model on two PULXs: M51 ULX-7 and NGC 7793 P13, using for both XMM-Newton observations. The aim was to test the validity of our model and derive information on the geometry of view, magnetic field, thermal and polarization properties of such sources.

In the model we considered the emission from a highly magnetized neutron star in an accreting binary system with a mass of 1.4 M_⊙, radius 10 km, and magnetic field strength between 10¹² and ∼10¹³ G. The accreting material originates from the donor star and proceeds through a geometrically-thin accretion disk up to the Alfvén radius, where particles are funneled along the magnetic field lines, which follow a dipolar field topology, towards the magnetic poles of the star. Each point of the torus and the disk emits like a blackbody, with a local temperature that depends on the magnetic colatitude θ, if the point belongs to the torus, or on the radial distance if the point belongs to the disk. The emission of the whole source is then a complex multicolor blackbody. We considered as well polarization observables, taking as polarized only the radiation coming from below the adiabatic radius r_pl, while that coming from the disk and the part of the torus outside the region of adiabatic propagation is unpolarized. We computed the flux and the polarization observables for different viewing geometries (χ, ξ) of each source.

We made direct comparisons of the light curves considering the 2–3 keV and 3–4 keV energy bands. We compared each observed profile with six simulations having different geometries of view, and from these comparisons, we constrained the viewing geometry obtaining a lower and upper limit of the angle χ, which was then used in the spectral analysis.

We compared spectra from different observations with those simulated with the model, varying the magnetic field strength, and using the temperature at the inner radius of the disk, T_in, as a free parameter. For M51 ULX-7 we used three XMM-Newton observations (0824450901, 0830191501 and 0830191601), and compared the spectrum with a multicomponent model that combines a blackbody (BBODY, between 0.1 − 1.5 keV), a multicolor blackbody (output of the model, between 1.5 − 5 keV), and a cut-off power law (CUTOFFPL, between 5 − 10 keV), assuming a magnetic field strength B = 8 × 10¹² G. The best fitting internal disk temperature is T_in ≈ 0.3 keV. We did the same for NGC 7793 P13, using the XMM-Newton observations 0693760401 and 0748390901, finding a T_in ≈ 0.32 keV.

For a given viewing geometry, the polarization degree increases with the energy and shows an oscillating behavior with the phase. The former is a consequence of the strong magnetic field considered in the model, while the oscillating behavior with the phase is a consequence of the rotation of the source combined with the geometry of view. Fixing the energy and phase but varying the angle χ we note that the value of the polarization degree increases. The variation of PD from one geometry of view to another could allow us to distinguish them if observed polarization measurements are available. Therefore, the polarization analysis has proved to be crucial in determining the geometry of view of the source when the spectral one gives degenerate results.

The model reproduces well the observed spectra of the two pulsating ULXs, particularly in the 1.5–5 keV interval, where the contribution of the torus with the disk is dominant.

Acknowledgments

SC is supported by the INAF Doctorate program. SC and LZ acknowledge financial support from the INAF Research Grant “Uncovering the optical beat of the fastest magnetized neutron stars (FANS)”. RTa and RTu acknowledge financial support from the Italian MUR through grant PRIN 2022LWPEXW. SC thanks E. Giunchi, B. Ambrosio and I. Salmaso for technical support.

References

All Tables

All Figures

	Fig. 1. A representation of the system. A 3D view on the top and a meridional cut on the bottom. The accretion disk and torus are shown in orange and blue respectively, the spin axis Ω and the magnetic dipole axis bdip are also indicated (in gray the radiation emitted from the accretion column).
In the text

	Fig. 2. LOS reference frame. The Z-axis is aligned to l, the LOS unit vector, η is the angle between the star magnetic axis bdip and the LOS, while ζ is the corresponding azimuth. The LOS and magnetic axis inclinations χ and ξ with respect to the spin axis Ω are also shown.
In the text

Fig. 3. The selection of the part in view for the torus, in panels (a) and (b), while considering also the disk, in panels (c) and (d). In panels (a) and (c) the source is shown exactly along the LOS direction, while in panels (b) and (d) the reference frame is rotated, for ease of visualization. The points in view are marked in cyan for the torus and in orange for the disk, while those not in view are marked in yellow for the torus and white for the disk. In this case, we used χ = 45°, ξ = 1°, Rmax = 70 RNS, and Rdisk = 200 RNS.

In the text

	Fig. 4. Comparison of the three XMM-Newton light curves of M51 ULX-7 with simulations. From top to bottom geometry of view with χ = 60° ξ = 10° (coral), with χ = 20° ξ = 10° (orange), and with χ = 30° ξ = 10° (green). The column on the left refers to the energy band 2–3 keV, while that on the right to the energy band 3–4 keV.
In the text

	Fig. 5. Comparison of the two XMM-Newton light curves of NGC 7793 P13 with simulations. From left to right the geometry of view with χ = 60° ξ = 10° (coral) and with χ = 40° ξ = 10° (red). The column on the right refers to the energy band 2–3 keV, while that on the left refers to the energy band 3–4 keV.
In the text

	Fig. 6. X-ray spectrum of M51 ULX-7 obtained with XMM-Newton (observation: 0824450901). The spectrum is fitted with a multicomponent model. The panel on the left refers to the viewing geometry χ = 20° ξ = 10°, that on the right to χ = 60° ξ = 10°. One component is phenomenological (BBODY), and one is calculated with the model (TORUSDISK).
In the text

	Fig. 7. X-ray spectra of M51 ULX-7 obtained with XMM-Newton (observations: 0824450901, 0830191501, 0830191601). The spectrum is fitted with a multicomponent model. The panels on the left refer to the viewing geometry χ = 20° ξ = 10°, while those on the right to χ = 60° ξ = 10°. Two components are phenomenological (BBODY, CUTOFFPL), and one is computed with the model (TORUSDISK).
In the text

	Fig. 8. X-ray spectrum of NGC 7793 P13 obtained with XMM-Newton (observations: 0693760401 on the top, 0748390901 on the bottom). The spectrum is fitted with a multicomponent model. Two components are phenomenological (BBODY, CUTOFFPL), and one is computed with the model (TORUSDISK). Panels on the top (bottom) from left to right are comparisons with geometries χ = 40° and χ = 60°, respectively.
In the text

	Fig. 9. PD and PA variation with energy and phase for different geometries of view. From left to right χ = 60°, 40°, and 20°.
In the text

	Fig. 10. PD and PA variation with energy and phase for different geometries of view. From left to right χ = 60° and 40°.
In the text

	Fig. 11. Polarization degree variation for different values of ΠL, for M51 ULX-7 with B = 8 × 1012 G and geometry of view χ = 60° ξ = 10°.
In the text

	Fig. 12. Polarization degree variation for different values of ΠL, for NGC 7793 P13 with B = 4 × 1012 G and geometry of view χ = 60° ξ = 10°.
In the text