Astronomy & Astrophysics manuscript no.

aanda ©ESO 2024

November 22, 2024

Pulsating ultraluminous X-ray sources: modeling the thermal

emission and polarization properties
S. Conforti1,2 , L. Zampieri2 , R. Taverna3 , R. Turolla3,4 , N. Brice6 , F. Pintore5 , and G.L. Israel7

1
Dipartimento di Fisica e Astronomia “G. Galilei”, Università degli Studi di Padova, Vicolo dell’Osservatorio 3, 35122 Padova, Italy
e-mail: [email protected]
2
INAF-Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35122 Padova, Italy
3
Dipartimento di Fisica e Astronomia “G.Galilei”, Università degli Studi di Padova, Via Marzolo 8, 35121 Padova, Italy
4
Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Surrey, RH5 6NT, UK
5
INAF/IASF Palermo, Via Ugo la Malfa 153, 90146 Palermo, Italy
arXiv:2411.13659v1 [astro-ph.HE] 20 Nov 2024

6
Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK
7
INAF–Osservatorio Astronomico di Roma, via Frascati 33, I-00078 Monte Porzio Catone, Italy
Received ...; accepted ...

ABSTRACT

Context. Ultraluminous X-ray sources (ULXs) are enigmatic sources first discovered in the 1980s in external galaxies. They are
characterized by their extraordinarily high X-ray luminosity, which often exceeds 1040 erg s−1 .
Aims. Our study aims to obtain more information about pulsating ULXs (PULXs), first of all, their viewing geometry, since it affects
almost all the observables, such as the flux, the pulsed fraction, the polarization degree (PD), and polarization angle (PA).
Methods. We present a simplified model, which primarily describes the thermal emission from an accreting, highly magnetized
neutron star, simulating the contributions of an accretion disk and an accretion envelope surrounding the star magnetosphere, both
described by a multicolor blackbody. Numerical calculations are used to determine the flux, PD, and PA of the emitted radiation,
considering various viewing geometries. The model predictions are then compared to the observed spectra of two PULXs, M51
ULX-7 and NGC 7793 P13.
Results. We identified the best fitting geometries for these sources, obtaining values of the pulsed fraction and the temperature at
the inner radius of the disk compatible with those obtained from previous works. We also found that measuring the polarization
observables can give considerable additional information on the source.
Key words. Neutron Stars – X-ray emission – Polarization

1. Introduction case, however, emission needs to be beamed to match the ob-

served luminosity (Fabrika & Mescheryakov 2001; King et al.
Ultraluminous X-ray sources (ULXs) were discovered in exter- 2001; Poutanen et al. 2007b).
nal galaxies in the 1980s with the Einstein X-ray Observatory
(Fabbiano 1989; Long et al. 1981). One of the most distinctive This picture was soon revolutionized in 2014, when pulsa-
features of ULXs is the high X-ray luminosity, which can reach tions were discovered in M82 X-2, a ULX in the M82 galaxy
values of L ∼ 1041 erg s−1 , far exceeding the Eddington limit for (Bachetti et al. 2014), supporting the view that at least some of
a solar mass object (see e.g. Kaaret et al. 2001, for a review). these sources could be powered by a neutron star (NS). The dis-
From an observational standpoint, ULXs are characterized by a covery of the first pulsating ULX was soon followed by several
higher occurrence in star-forming galaxies (Swartz et al. 2004), others (Fürst et al. 2016; Israel et al. 2017a,b; Rodríguez Castillo
a turnover in the spectra at energies in the 2–10 keV range, a et al. 2020; Carpano et al. 2018; Sathyaprakash et al. 2019;
second softer component below 1 keV, and spectral variability Quintin et al. 2021), highlighting the existence of an entire pop-
between different epochs (e.g. Gladstone et al. 2009; Pintore & ulation of pulsating ULXs, dubbed PULXs.
Zampieri 2012; Sutton et al. 2013a,b; Middleton et al. 2015a,b). Even though the association of PULXs with neutron stars
The large amount of energy emitted and the X-ray variabil- appears well established, explaining how neutron stars can emit
ity led to believe that ULXs were binary systems in which a such a huge amount of energy remains an open problem, given
massive donor transfers material onto a black hole (Kaaret et al. that the Eddington luminosity for a neutron star is ≈ 1038 erg s−1 .
2017). For some time, ULXs were thought to be powered by Accretion onto (highly) magnetized NSs has been the focus of
black holes of intermediate mass (IMBHs), that is, larger than many investigations, starting from the pioneering work of Basko
100M⊙ (Colbert & Mushotzky 1999; Makishima et al. 2000; & Sunyaev (1976). They presented a model for accretion onto
Miller et al. 2003). In fact, the Eddington luminosity for such magnetized neutron stars that included the formation of an ac-
objects exceeds the observed one, so accretion can proceed at cretion column on the magnetic poles of the star for the first
a sub-Eddington rate. Later on, super-Eddington accretion onto time. This model for the column was somewhat simplistic and
stellar-mass black holes (MBH ∼ 10–20 M⊙ , Stobbart et al. 2006; did not take into account the consequences of the propeller ef-
Feng & Soria 2011 and also black holes with MBH ∼ 30–80 M⊙ , fect, induced by the huge magnetic field strength (∼ 1015 G),
Zampieri & Roberts 2009) were also considered. If this was the on the column. Subsequent investigations addressed the prob-
Article number, page 1 of 17
 A&A proofs: manuscript no. aanda

lem with an increasing level of sophistication, including the ef- we provide a brief description of how the data were analyzed
fects of a super-strong magnetic field (see Lyubarskii & Syun- while section 4 is devoted to the spectral and polarization results,
yaev 1988; Mushtukov et al. 2015, 2017; Brice et al. 2021). In which are then discussed in section 5. Finally, our conclusions
particular, it was realized that in a strongly magnetized accre- are summarized in section 6.
tion 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, al- 2. Physical model
lowing for a larger flux to escape from the sides of the column. In this section, we introduce the physical model adopted to sim-
Such a scenario could indeed explain the large luminosity ob- ulate and reproduce the X-ray spectrum and polarization proper-
served from these sources. We present here a simplified model ties of the radiation emitted by pulsating ULXs.
(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 2.1. Modeling the source
by the magnetic field lines that extend up to the magnetospheric The model computes the emission from an accreting magnetized
radius. To this aim, we adapted the numerical code by Taverna neutron star (NS) in a binary system. The accreting material orig-
& Turolla (2017), which calculates the flux of photons coming inates from the donor star and proceeds through a geometrically-
from the part in view of the system as a function of both the pho- thin accretion disk. At the Alfvén radius (rA ), where the mag-
ton energy and the star rotational phase, so to compute spectra netic pressure equals the ram pressure of the gas particles, the
and light curves of PULXs in different energy bands. In addition matter is funneled along the closed field lines of the magnetic
to the flux, the code also computes the polarization degree (PD) field, falling towards the polar regions. This envelopes the star in
and polarization angle (PA) as functions of energy and phase. a sort of cocoon that, for the sake of simplicity, we call “torus”
Finally, we compare the simulations with the observational hereafter (Mushtukov et al. 2015, 2017; Brice et al. 2021). At
data of two PULXs, to test if the model can reproduce the prop- the high accretion rates expected to occur in PULXs, this torus
erties exhibited by real sources. Firstly we focus on M51 ULX- is optically thick. In the following, we assume that the neutron
7 (aka NGC 5194 X-7, Roberts & Warwick 2000, CXOM51 star has a mass of 1.4 M⊙ , a radius of 10 km, and that the mag-
J133001.0+47134, Terashima & Wilson 2004, NGC 5194/5 netic field has a dipolar topology with a strength (at the pole) in
ULX-7, Liu & Mirabel 2005), located in the outskirts of a young the range 1012 –1013 G. The equation of the magnetic field lines
open cluster in a spiral arm of its host galaxy. It belongs to a is, in polar coordinates,
high-mass X-ray binary (HMXB) (Abolmasov et al. 2007) with
a (likely) O-B companion with mass ≳ 8 M⊙ (Rodríguez Castillo r = rA sin2 θ , (1)
et al. 2020; Abolmasov et al. 2007) and exhibits an X-ray lumi-
nosity of 6 × 1039 erg s−1 . It was suggested that the source is where θ is the magnetic colatitude, and rA , for the sake of sim-
powered by a neutron star accreting at a super-Eddington rate plicity, here coincides with the maximum radius of the magne-
(Brightman et al. 2022). NGC 7793 P13 (P13 hereafter) is a pul- tosphere, which is computed with the equation (1) in Brice et al.
sating ULX located on the southern edge of the galaxy NGC (2023):
7793 (in the Sculptor group Fabbiano et al. 1992), the X-ray Rm = 7 × 107 ΛM 1/7 R10/7 B4/7 −2/7
6 d,12 L39 cm, (2)
emission of which is completely dominated by this source. The
ESO-VLT observations in 2008–2009 revealed that NGC 7793 where Λ is a dimensionless parameter that depends on the mode
P13 has a B9Ia supergiant companion with a mass between 10 of accretion (a typical value of 0.5 for accretion via thin disc),
and 20 M⊙ (Motch et al. 2011; Fürst et al. 2021). From XMM- Bd,12 is the surface dipole field strength at the magnetic poles
Newton and NuSTAR observations in 2016, coherent pulsations (expressed in units of 1012 G), and L39 is the accretion lumi-
were detected (P ≈ 0.42 s); this led to identify P13 as the third nosity (in units of 1039 erg s−1 ). The equation (1) is just an ap-
known ULX powered by a neutron star (Israel et al. 2017b; Fürst proximation since the exact topology of the magnetic field near
et al. 2016). the surface of the NS is unknown. Pulsating ULXs, in particular,
The main goal of this paper is to derive information on the may have a multipolar magnetic field structure, as proposed by
viewing geometry of these sources confronting the simulated Israel et al. (2017a) to accommodate for the large field required
light curves with the observed ones, to constrain the temperature to rise the Eddington limit enough and, at the same time, prevent
at the inner radius of the disk, and the magnetic field strength. the star from entering the propeller stage. However, higher-order
To reproduce the spectrum in the 0.1 − 10 keV band, three com- multipoles quickly decay moving away from the surface, so our
ponents are required: the spectrum from our model (torus and assumption is quite justified if rA is tens stellar radii and makes
disk), a cut-off power law component to account for the non- the computations much simpler.
thermal emission in the hard band (believed to come from the Both the disk and the torus are assumed to be in Local Ther-
accretion column), and a soft blackbody component, possibly mal Equilibrium (LTE). The disk temperature follows the thin-
due to radiative winds produced in the intermediate/outer part of disc profile (see Shakura & Sunyaev 1973)
the accretion disk (Poutanen et al. 2007b; Pinto et al. 2016, 2020;  R 3/4
Middleton et al. 2022). In case spectral fitting is degenerate for T (r) = T in
in
, (3)
different geometries, we use polarization to discriminate among r
them. This further shows that polarimetric measurements are in- where T in is the temperature at the inner radius of the disk, Rin ,
deed a powerful tool for the study of highly magnetized neutron that we take to coincide with the radius at which the torus and
stars in PULXs, as testified by the results obtained by the recent the disk intersect. The temperature profile of the torus is that
NASA-ASI mission IXPE (Imaging X-ray Polarimetry Explorer, presented in Brice et al. (2023), it varies with the magnetic co-
Weisskopf et al. 2022). latitude and follows the expression:
This work is organized as follows. In section 2 we present
the model and the assumptions on which it relies. In section 3 T out,torus = T in,torus τ−1/4 , (4)
Article number, page 2 of 17
 S. Conforti et al.: Pulsating ultraluminous X-ray sources: modeling the thermal emission and polarization properties

k and the local magnetic field B, while in the latter the elec-
tric field oscillates perpendicularly to both these vectors. The
cross-sections of X-mode photons are greatly reduced below the
electron cyclotron energy, Ec,e ≈ 11.6(B/1012 G) keV, by a fac-
tor of ∼ (B/BQ )2 with respect to the unmagnetized case, where
BQ ≃ 4.414 × 1013 G is the critical field (Mészáros 1992; Hard-
ing & 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 mod-
ify the optical properties of the medium in which radiation prop-
agates. In particular, vacuum birefringence (Heisenberg & Euler
1936) affects photons propagating in a strongly magnetized vac-
uum, 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 in-
stantaneously 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) ra-
dius (Heyl et al. 2003; Taverna et al. 2015):
 E 1/5  R 1/5 !2/5
rpl NS Bp
Fig. 1: A representation of the system. A 3D view on the top and ∼5 , (5)
RNS 1 keV 10 km 1011 G
a meridional cut on the bottom. The accretion disk and torus are
shown in orange and blue respectively, the spin axis Ω and the where E is the photon energy, Bp the polar magnetic field
magnetic dipole axis bdip are also indicated (in gray the radiation strength, and RNS the NS radius. The actual polarization state
emitted from the accretion column). 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 < rpl while
where T out,torus is the temperature at the torus surface, T in,torus the the electric field is frozen for r > rpl . For the magnetic field
temperature at the torus inner boundary, and τ > 1 the torus strength (∼ 1013 G) and the photon energies considered here
optical depth (see figure 1 in Brice et al. 2023). As reported in (0.1–1 keV), the Alfvén radius (rA ) is larger than the adiabatic
Brice et al. (2023), the local optical depth depends on the dy- radius (rA ∼ 50RNS , while rpl ∼ 36RNS ), so we expect that a
namics of the accreting particles: the non-uniform velocity field large part of the torus lies outside the adiabatic region. Hence,
inside the envelope introduces a meridional variation of T out . Be- we assume radiation coming from this zone of the torus to be un-
cause of our assumption of LTE, each annulus of both the disc polarized1 . Radiation coming from inside rpl carries a non-zero
and the torus emits blackbody radiation at the local temperature; intrinsic polarization degree,
no effects due to the reprocessing of thermal radiation by mate-
rial on top of the disc and/or torus are considered. In our sim- NX − NO
ΠL = , (6)
plified modeling, we do not account for radiation coming from NX + NO
the sides of the accretion columns which form above the mag-
netic poles if the mass accretion rate is high enough. Emission where NX (NO ) is the fraction of X (O) photons. Following the
from the columns gives rise to a non-thermal tail in the X-ray same argument, we assume that the polarization of radiation
spectrum (Walton (2015),Walton et al. (2018a), Brightman et al. coming from the disk is negligible since the disk starts outside
(2022)). In order to reproduce the observed spectrum over the rpl where B is quite small2 .
entire range, we added to the model a phenomenological cut-off It is convenient to express the polarization observables, the
power law, often used to describe this type of non-thermal emis- polarization degree (PD) and polarization angle (PA), in terms of
sion (see e.g. Walton et al. 2018b) the Stokes parameters I, Q, and U (Rybicki & Lightman 1991).
In a reference frame with the z axis along the unit wavevector k,
2.2. Polarization 1
This is justified because, whatever the intrinsic polarization is, the
polarization degree measured at infinity is close to zero as a conse-
Radiation emitted by strongly magnetized neutron stars is ex- quence of geometric effects (Taverna et al. 2015, see also below).
pected to be highly polarized in two normal modes, the ordinary 2
We are interested in the radiation coming from the torus that is polar-
(O) and extraordinary (X) ones (e.g. Gnedin & Pavlov 1974; Ho ized in two normal modes, that depend on the magnetic field direction.
& Lai 2003; Lai et al. 2010). The former is characterized by an So, for the sake of simplicity, we neglect the polarization of the disk
electric field oscillating in the plane of the propagation vector radiation, since it is different from that of the torus.

Article number, page 3 of 17

 A&A proofs: manuscript no. aanda

and the y axis in the (k, B) plane, the Stokes parameters for X
and O photons are
 I  1  I   1 
       
 Q  = 1  Q  = −1 . (7)
U X 0 U O 0
The last Stokes parameter, V, which describes the circular po-
larization, is not considered here. The Stokes parameters de-
pend on the reference frame of each photon (xi , yi , zi ), which is
fixed by the local direction of the (dipole) magnetic field that
changes across the emitting surface. To sum the Stokes param-
eter 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 zi ) with respect to the polarimeter frame
by an angle αi . Eventually, for a single photon the Stokes param-
eters Ii , Qi , Ui are defined as
Ii = Ii
Qi = Qi cos(2αi ) + Ui sin(2αi ) (8) Fig. 2: LOS reference frame. The Z-axis is aligned to l, the LOS
Ui = Ui cos(2αi ) − Qi sin(2αi ). unit vector, η is the angle between the star magnetic axis bdip
The observed PD and PA are given by (Rybicki & Lightman and the LOS, while ζ is the corresponding azimuth. The LOS
1991) and magnetic axis inclinations χ and ξ with respect to the spin
axis Ω are also shown.
Q2 + U 2
p
PD =
I (9)
2.3.1. Visibility of the source
!
1 U
PA = arctan ,
2 Q The code starts by computing the visible part of the source as-
where the Stokes parameters Q and U refer to the “total” radia- suming a specific viewing geometry. We consider the surface,
tion, that is, the sum of the Qi and Ui of the collected photons. comprising the disk and the torus, as a collection of small emit-
Substituting equations (7) and (8) inside equations (9) one ob- ting patches. The grid is made by 50 bins in θ, ϕ, and the radial
tains distance r, while for both the energy and rotational phase we take
1h XX 30 bins. The phase is measured from the projection of the line of
PD = N+2 cos (2αi − 2αk ) sight (LOS) onto the plane perpendicular to the spin axis.
N i k>1 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
XX
+2 cos (2α j − 2αr )
j r> j for the chosen angles χ and ξ, which represent the orientations
XX i1/2 (10) of the neutron star spin axis and the magnetic axis with respect
−2 cos (2αi − 2α j ) to the LOS, respectively. The code also takes into account the
i j orientation of the disk relative to the NS spin axis. The follow-
ing steps are performed: selecting the visible part of the torus,
P P
1 i sin(2αi ) − j sin(2α j )
PA = arctan − P P , accounting for its self-shadowing as well as the shadow cast by
2 i cos(2αi ) − j cos(2α j )
the disk, and selecting the visible part of the disk not shadowed
where i, k = 1, ..., NX and r, j = 1, ..., NO . From equations (10) by the torus.
one can see that the polarization degree measured at infinity is To achieve this, the code traces a straight line parallel to the
in general lower than that at the emission (equation 6), because LOS and identifies the intersections between this line and the
of the presence of the geometrical factors sin(2α) and cos(2α). If emitting surface (torus and/or disk). To understand if each inter-
the magnetic field topology is particularly tangled (as it occurs section is in view, we evaluate the quantity
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 n· l > 0, (11)
field (like that at rpl , far from the surface), the different α angles
attain much more similar values and radiation is much less depo- where n represents the surface normal unit vector at the inter-
larized (Taverna et al. 2015). Hence, when vacuum birefringence section, and l is the unit vector of the LOS. If multiple intersec-
is accounted for and the Stokes parameters are computed at rpl , tions satisfy equation (11), we calculate the z-coordinate of each
the polarization pattern at the surface is likely preserved also at point and consider in view the one with the largest value of the
the observer. z-coordinate.
Figure 3 shows the emitting regions of the source for a given
2.3. Numerical implementation viewing geometry, in the LOS reference frame. Although a full
treatment requires to take into account general relativistic (GR)
To compute the spectral and polarization properties of the emit- effects, such as gravitational redshift and ray-bending, we expect
ted radiation, we used the ray-tracer code discussed in Taverna them to be negligible since the emitting regions (torus and disc)
et al. (2015), but adapted to our model, in particular to a different lie in general far from the star. For these reasons, we do not con-
emitting surface (torus and disk). sider GR corrections in the present simulations.
Article number, page 4 of 17
 S. Conforti et al.: Pulsating ultraluminous X-ray sources: modeling the thermal emission and polarization properties

(a) (b)

(c) (d)

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 tion 3.2), obtaining as output a matrix with the total observed
flux as a function of energy and phase.
Once a point within the emission region of the torus is tagged
as in view, we calculate its temperature. As reported above, the 2.3.3. Polarization observables
temperature of the torus varies with the magnetic colatitude. The
temperature of each emission point is obtained by interpolat- For each point in view inside rpl we set NX (NO = 1 − NX )
ing over the grid discussed in Brice et al. (2023). Instead, for and compute the rotated Stokes parameters of the photons at
a visible point on the disk we derive the temperature from the the adiabatic radius (equation 5 in Taverna et al. 2015). As de-
Shakura-Sunyaev profile (equation 3), once the inner radius Rin scribed above for the total flux, we sum the Stokes parameter
and the inner temperature T in of the disk are fixed (subsection fluxes over the part in view of the emitting region (at each en-
4.2). The locally emitted flux is calculated from the blackbody ergy and phase), obtaining the total Q and U fluxes for both X
distribution at the local temperature. We then sum all the con- and O photons. Hence, the code calculates PD and PA following
tributions coming from different points in view using the same equation (9). For the points in view that are at a distance ≳ rpl
approach as discussed in Taverna & Turolla (2017, see their sec- from the star we do essentially the same. However, since, as de-
Article number, page 5 of 17
 A&A proofs: manuscript no. aanda

scribed above, outside the adiabatic region the photon electric simulations having different geometries of view, each differing
field direction is frozen with respect to the magnetic field one, only for the angle χ (ranging from 10◦ to 60◦ ). The value of ξ was
we anyway expect to observe a much lower polarization degree fixed at 10◦ , because for values of 20◦ , or higher, and for values
for radiation emitted from these points. For this reason, and for lower than 10◦ , the PFs were not compatible with that observed.
the sake of simplicity, for photons coming from outside rpl we So, for the sake of simplicity, we fixed ξ at 10◦ , also to represent
fixed directly NX = 0.5 (considering the 50% of radiation to be a general non-aligned magnetic rotator case. χ angles beyond
X-mode photons, and the other 50% to be O-mode photons, with 60◦ were not considered, because we noted in the simulations
no privileged polarization mode), which, from equation (6) gives that the PF decreases if χ > 60◦ .
ΠL = 0. So, in this work, we analyzed the variation of PA and
PD with the energy and phase varying NX only for photons com-
ing from the regions of the torus closer to the star surface (and 4.1.1. M51 ULX-7
so emitted from inside rpl ). To perform simulations we fixed some input parameters: the
magnetic field strength B, the Alfvén radius rA , the torus tem-
perature, and the disk radius Rdisk . M51 ULX-7 has a surface
3. Data reduction magnetic field strength between 8 × 1011 G and 1013 G (values
We analyzed a sample of XMM-Newton and NuSTAR observa- derived by Rodríguez Castillo et al. 2020, from the P-Ṗ relation),
tions of M51 ULX-7 and NGC 7793 P13, two PULXs taken as so we set a lower and upper limit, B = 1012 G and B = 8 × 1012
benchmarks for our study. We used three XMM-Newton obser- G, respectively. The Alfvén radius (rA , which coincides with the
vations for M51 ULX-7 (Obs.ID: 0824450901, 0830191501 and inner radius of the disk Rin ) and the temperature of the torus were
0830191601), carried out between 2018 May and June, and two calculated according to B. The outer radius of the disk was set
for NGC 7793 P13 (Obs.ID: 0693760401 and 0748390901), car- to 200RNS , as beyond this radius the emitted flux falls below 0.1
ried out between 2013 November and 2014 December. We ex- keV, outside the range considered here. We report all the input
tracted the EPIC-pn data using the SAS v.14.0. package3 . Spec- fixed parameters in table 1.
tra were obtained by selecting events with the pattern ≤ 4 for For M51 ULX-7, we found a lower and upper limit for the
EPIC-pn (single- and double-pixel events) and setting ‘flag=0’ geometry of view, χ = 20◦ ξ = 10◦ and χ = 60◦ ξ = 10◦ ,
to ignore bad pixels and events coming from the CCD edges. respectively (see figure 4). In table 2 we report the values of
Epochs of high background were also removed from the analy- the observed PFs with those of the best-fit viewing geometries,
sis. Source and background events were extracted from circular across the two different energy bands. The observed PFs show
regions with radii of 30′′ and 60′′ , respectively. The spectra were general agreement with the simulated ones within uncertainties.
rebinned to have at least 25 counts per energy bin. In figure 4 we show the comparisons.

B (G) T torus,max (keV) T torus,min (keV) rA (RNS )

4. Results
12
10 1.38 1.06 19
We started by comparing the light curves of the XMM-Newton
observations with the simulated ones to constrain the geome- 8 × 1012 0.77 0.57 52
try of view, which significantly affects the flux modulation. In
particular, we compared the observed and simulated pulsed frac- Table 1: Magnetic filed strength B, minimum and maximum tem-
tions (PFs), defined as perature of the torus, T torus,min and T torus,max respectively, and the
Alfvén radius, rA , for spectral simulations of M51 ULX-7.
max(Flux) − min(Flux)
PF = . (12)
max(Flux) + min(Flux)
After constraining the viewing geometry, we compared the spec- 4.1.2. NGC 7793 P13
tra determining the temperature at the inner radius of the accre-
tion disc (T in ), the torus temperature range, and the magnetic We did the same for NGC 7793 P13, considering that the ob-
field strength. For the comparison of the pulse profiles, we con- served surface magnetic dipole field strength is around 2 × 1012
sidered six possible viewing geometries while, for the compar- G (Fürst et al. 2016; Israel et al. 2017b, still obtained from
ison of the spectra, we chose 16 values of T in , selected on the P and Ṗ). We then set the dipole magnetic field strength to
basis of values found in the literature, for each different value B = 4 × 1012 G considering the value at the poles, and computed
of the magnetic field strength. We then analyzed the polarization the Alfvén radius of the torus, the torus temperatures (see table
degree and polarization angle for each source using the best-fit 3), and set the disk outer radius to Rdisk = 200 RNS . This time we
viewing geometry, T in , and magnetic field strength. found much less variation for the angle χ, identifying χ = 60◦
ξ = 10◦ and χ = 40◦ ξ = 10◦ as the upper and lower limit, re-
spectively. In figure 5 the light curves and in table 4 the pulsed
4.1. Comparison of the light curves fractions, across the two energy bands, are confronted with the
We made a direct comparison of the light curves, considering the observed ones.
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 4.2. Comparison of the spectra
torus) dominates. We compared each observed profile with six
4.2.1. M51 ULX-7
3
Spectra reduced with this version were already available in our
archives. Since they are of enough high quality for the aims of this work In section 4.1.1 we derived a lower and upper limit for the ge-
and since we do not expect significant differences, we choose not to re- ometry of view of the source, χ = 20◦ , ξ = 10◦ and χ = 60◦ ,
process the data with the latest SAS v21.0.0. ξ = 10◦ , respectively. We now perform spectral comparisons
Article number, page 6 of 17
 S. Conforti et al.: Pulsating ultraluminous X-ray sources: modeling the thermal emission and polarization properties

Obs. ID Observed PFs (%) Simulated PFs (%)

2–3 keV 3–4 keV 2–3 keV 3–4 keV
0824450901 2631
22 3227
37 26.48d 29.2d
0830191501 1216
8 1217
7 14.2a 15.0a
0830191601 1721
13 1823
13 20.1b 21.5b

Table 2: Observed pulsed fractions of M51 ULX-7 for the three XMM-Newton observations, against simulated ones both calculated
in the energy bands 2–3 keV and 3–4 keV.
Notes.
(a)
Geometry of view with χ = 20◦
(b)
Geometry of view with χ = 30◦
(c)
Geometry of view with χ = 60◦

B (G) T torus,max (keV) T torus,min (keV) rA (RNS ) 4.2.3. M51 ULX-7 – B = 8 × 1012 G

4 × 1012 0.79 0.19 49 In this case, we reproduced the spectrum by the superposition of
a soft X-ray thermal component plus a hard X-ray non-thermal
Table 3: Magnetic filed strength B, minimum and maximum tem- one. We fitted the former with the sum of two thermal spec-
perature of the torus, T torus,min and T torus,max , respectively, and the tral components, the XSPEC bbody model and the torusdisk
Alfvén radius, rA , for spectral simulations of NGC 7793 P13. 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 ana-
to constrain T in for both the geometries and the magnetic field lytical expressions of cutoffpl. Still to ensure the convergence of
strengths. From the literature, the temperature of the soft spec- the fit, we fixed the temperature of the bbody and also the column
tral component (which is related to T in ) varies between 0.33 and density, NH = 9 × 1019 cm−2 5 .
0.5 keV (Rodríguez Castillo et al. 2020; Brightman et al. 2022; We used the superposition of a blackbody, a multicolor
Vasilopoulos et al. 2020). We then produced 16 spectral simula- blackbody and a power law component as done by Brightman
tions, each with a different temperature value T in , which ranges et al. (2022). The low-energy blackbody component, in fact, is
between 0.25 and 0.5 keV. used to describe the emission due to radiative winds that can
originate from the accretion disk at super-Eddington accretion
4.2.2. M51 ULX-7 – B = 1012 G rates, (see e.g. Poutanen et al. 2007a; Pinto et al. 2016, 2020).
In figure 7 we show the spectra for the three XMM-Newton ob-
We reproduced the spectrum with the superposition of two ther- servations. In table 6 we report the best-fit values of the parame-
mal components, the XSPEC bbody model (see Arnaud 1996) ters. The value of T in turns out to be ∼ 0.3 keV, compatible with
and our simulated multicolor blackbody, similarly as in Ro- those on the literature, for both the viewing geometries, which
dríguez Castillo et al. (2020). We considered the interstellar also show indistinguishable spectral results. Considering the re-
absorption through the wabs model (Morrison & McCammon sults with B = 8 × 1012 G, in the following, we took into account
1983) leaving the column density as a free parameter of the fit. only the case with this value of the magnetic field strength.
We carried out the fitting procedure outside XSPEC, using a
Python script where we implemented the analytical expressions 4.2.4. NGC 7793 P13
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 We adopted the same approach for the spectral fits on NGC 7793
wabs and not tbabs, as an analytical expression is available for P13, using the geometries χ = 60◦ ξ = 10◦ and χ = 40◦ ξ = 10◦
the former. We fixed the emission temperature of the bbody at derived form the pulse profiles analysis. From the X-ray spectral
0.25 keV, to ensure the convergence of the fit4 . Furthermore, we analyses of Walton et al. (2018a), the spectrum of this source
stress that our purpose is to test the reliability of the torusdisk can be fitted with a model made by three components, two of
model by varying the contribution of the internal disk temper- which are thermal and dominate below 10 keV. They found that
ature T in for a given viewing geometry, while accounting for the temperature of the coolest thermal component spans from
the presence of additional components. In figure 6 we show the ∼ 0.3 to 0.5 keV, and that of the hottest thermal component from
spectral fits for both the geometries of view of the first obser- ∼ 1 to 1.5 keV. We used the same multicomponent model de-
vation (the others provide similar results). In table 5 we report scribed above (bbody + torusdisk model + cutoff power law),
the best-fit parameters of all the observations. The best value for correcting for absorption (wabs) with a fixed column density,
the temperature T in turns out to be 0.5 keV, in agreement with NH =∼ 3 × 1020 cm−2 , and temperature of the bbody compo-
Rodríguez Castillo et al. (2020) within the uncertainties. Despite nent, 0.25 keV. We then varied T in between 0.25– 0.5 keV in 16
the superposition of these two thermal components can fit the bins. We show the results in figure 8 and we report the values
data, it also leaves non negligible residuals between 1.5–3 keV. 5
To obtain this value of the column density, NH , we fitted singularly
the three spectral components, each in the energy range where they
4
The value of the bbody temperature is calculated by separately fitting dominate. Then we calculated NH as the mean of the three values ob-
this component in the energy range where it dominates. tained from the fits.

Article number, page 7 of 17

 A&A proofs: manuscript no. aanda

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.

of the parameters in table 7. We find T in ∼ 0.3 keV, in agree- ing from the adiabatic region, we fixed the intrinsic polarization
ment with the literature, and also in this case the spectra of the fraction, ΠL , at 60% in the X-mode. We did not computed this
two geometries of view do not show differences. The results also value self-consistently since this would have required to solve
suggest that B = 4 × 1012 G can reproduce the spectrum well. the full radiation transport problem, which is definitely outside
the scope of the present work. Nevertheless, in the presence of
strong magnetic fields (≳ 1013 G) X-mode photons are expected
4.3. Polarization to dominate the emission, so our assumption is not unreasonable.
The effects of changing ΠL will be discussed further on. We then
In simulating the polarization properties of the two PULXs, we
analyzed the variation of the polarization degree and polariza-
considered fields with strengths of 4 × 1012 G (NGC 7793 P13)
tion angle with the geometry of view. In the case of M51 ULX-7
and 8 × 1012 G (M51 ULX-7) at the poles. For photons escap-
Article number, page 8 of 17
 S. Conforti et al.: Pulsating ultraluminous X-ray sources: modeling the thermal emission and polarization properties

Obs. ID Observed PFs (%) Simulated PFs (%)

2–3 (keV) 3–4 (keV) 2–3 (keV) 3–4 (keV)
0693760401 2326
19 3237
28 24.3a 29.8b
0748390901 2225
20 2932
26 23.7a 29.9b

Table 4: Observed pulsed fractions of NGC 7793 P13 for the two XMM-Newton observations against simulated ones both calculated
in the energy bands 2–3 keV and 3–4 keV.
Notes.
(a)
Geometry of view with χ = 40◦
(b)
Geometry of view with χ = 60◦

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.

we used χ = 60◦ , 40◦ and 20◦ , while for NGC 7793 P13 χ = 60◦ Fixing energy and phase but varying the angle χ from 20◦
and 40◦ . to 60◦ the value of the polarization degree increases. This effect
could be caused by the geometrical depolarization. The polar-
ization angle oscillates with the phase remaining almost constant
4.3.1. M51 ULX-7 with the energy. As for the polarization degree, the variation with
the phase of the polarization angle depends on the viewing ge-
We show in figure 9 the energy- and phase-dependent PD, for ometry.
three different viewing geometries. In each panel, the polariza-
tion degree increases with the energy and shows an oscillating
behavior with the phase, that changes going from χ = 20◦ to 4.3.2. NGC 7793 P13
χ = 60◦ . The former is a consequence of the strong magnetic
field considered in the model, while the latter is a consequence In figure 10 we show the polarization results for NGC 7793 P13.
of the rotation of the source combined with the geometry of view. The maximum polarization degree is lower than that for M51
Article number, page 9 of 17
 A&A proofs: manuscript no. aanda

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).

Table 5: Best fit parameters for M51 ULX-7 (B = 1012 G)

χ = 20◦ , ξ = 10◦
Observation ID a
N (10 ) 23
J (10−3 )
b
T in c NH χ2 /dof
(keV/cm2 /s) (keV/cm /s) (keV) (1022 cm−2 )
2

0824450901 4.17 ± 0.07 4.8 ± 0.2 0.5 3.3 ± 1.4 629/386

0830191501 3.17 ± 0.06 6.6 ± 0.3 0.5 4.0 ± 1.6 647/382
0830191601 3.27 ± 0.06 4.5 ± 0.2 0.5 2.4 ± 1.4 619/382

χ = 60◦ , ξ = 10◦
Observation ID Na (1023 ) Jb (10−3 ) T in c NH χ2 /dof
(keV/cm2 /s) (keV/cm /s) (keV) (10 cm−2 )
2 22

0824450901 5.73±0.09 4.7 ± 0.2 0.5 2.8 ± 1.4 558/386

0830191501 4.37 ± 0.08 6.5 ± 0.3 0.5 3.5 ± 1.5 584/382
0830191601 4.51 ± 0.07 4.4 ± 0.2 0.5 1.8 ± 1.3 555/382

Notes.
(a)
torusdisk norm
(b)
bbody norm
(c)
Temperature at the inner radius of the accretion disk

ULX-7. Fixing energy and phase, it is higher for the geometry the hypothesis that, even if the intrinsic polarization value as-
of view with χ = 40◦ than that with χ = 60◦ . As for M51 ULX- sumed in the model may be different from the one computed
7, the polarization degree oscillates with the phase and increases self-consistently, we can still distinguish the geometry of view
with the energy, while the polarization angle oscillates with the through the energy and phase variations of PD, since they do not
phase remaining quite constant with the energy for low values of depend on the value of ΠL , but only on the geometry of view and
χ. magnetic field of the source.

4.3.3. Polarization results with different ΠL 5. Discussion

We previously mentioned that the value of ΠL is not computed, This paper aims to constrain the geometry of view, the mag-
but it is fixed considering the magnetic field strength of the netic field strength, and thermal properties of pulsating ULXs,
source. We report here the results for the polarization degree reproducing the thermal radiation emitted by a highly magne-
varying ΠL and fixing the geometry of view with χ = 60◦ and tized neutron star accreting at a super-Eddington rate and taking
ξ = 10◦ , for both M51 ULX-7 and NGC 7793 P13. The re- into account the emission from an accretion disk and an opti-
sults show that a variation in the intrinsic polarization degree cally thick envelope surrounding the magnetosphere. The model
only leads to an overall variation in the value of the PD, with- reproduces the light curves and the spectra of PULXs, for the lat-
out changing its behavior with energy and phase. This reinforces ter by adding one or two phenomenological components that are
Article number, page 10 of 17
 S. Conforti et al.: Pulsating ultraluminous X-ray sources: modeling the thermal emission and polarization properties

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).

not self-consistently included in the model: a blackbody compo- 5.1. M51 ULX-7
nent to model the excess of flux at low energies, probably pro-
duced by optically thick winds from an intermediate region of 5.1.1. Light curves and spectra
the disk; a cut-off power law produced by an accretion column.
In addition to the spectral properties, we also incorporated the We started comparing the light curves because the variation of
polarization observables of the emitted radiation, focusing our the flux with the phase, and the corresponding pulsed fraction,
analyses on the polarization degree. We discuss our results in are highly dependent on the viewing geometry, allowing us to
the following. 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 ac-
cretion 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 solu-
tion for the viewing geometry but we were able to derive a lower
Article number, page 11 of 17
 A&A proofs: manuscript no. aanda

Table 6: Best fit parameters for M51 ULX-7 (B = 8 × 1012 G)

χ = 20◦ , ξ = 10◦
Observation ID a
N (10 ) 23 b
K (10 )−6
αb βb Jc (10−3 ) T in d χ2 /dof
(keV/cm2 /s) (keV/cm2 /s) (keV) (keV/cm2 /s) (keV)

0824450901 5.8 ± 0.8 1.73 ± 1.73 -5.1 ± 1.1 1.3 ± 0.3 3.3 ± 0.2 0.32 426/384
0830191501 3.7 ± 0.2 0.094 ± 0.093 -8.2 ± 1.1 0.9 ± 0.1 4.2 ± 0.2 0.32 445/380
0830191601 9.8 ± 7.6 39.3 ± 31.9 -1.6 ± 0.6 5.1 ± 1.9 3.6 ± 1.0 0.28 423/380

χ = 60◦ , ξ = 10◦
Observation ID Na (1023 ) Kb (10−6 ) αb βb Jc (10−3 ) T in d χ2 /dof
(keV/cm2 /s) (keV/cm2 /s) (keV) (keV/cm2 /s) (keV)

0824450901 10.1 ± 1.2 1.26 ± 1.28 -5.3 ± 1.1 1.3 ± 0.3 3.5 ± 0.2 0.32 431/384
0830191501 6.6 ± 0.3 0.046 ± 0.049 -8.9 ± 1.2 0.8 ± 0.1 4.5 ± 0.2 0.32 448/380
0830191601 19.2 ± 13.2 43.5 ± 26.9 -1.5 ± 0.4 5.6 ± 1.6 3.5 ± 0.9 0.30 420/380

Notes.
(a)
torusdisk norm
(b)
cutoffpl: norm, photon index, and e-folding energy of exponential cut-off
(c)
bbody norm
(d)
Temperature at the inner radius of the accretion disk

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.

Article number, page 12 of 17

 S. Conforti et al.: Pulsating ultraluminous X-ray sources: modeling the thermal emission and polarization properties

Table 7: Best fit parameters for NGC 7793 P13

χ = 40◦ , ξ = 10◦
Observation ID a
N (10 )23 b
K (10 ) −6
α b
βb Jc (10−3 ) T in d χ2 /dof
(keV/cm2 /s) (keV/cm2 /s) (keV) (keV/cm2 /s) (keV)

0693760401 3.5 ± 0.2 0.6 ± 0.4 -7.5 ± 0.8 0.8 ± 0.1 4.5 ± 0.3 0.32 329/197
0748390901 1.46 ± 0.08 2.0 ± 1.1 -6.5 ± 0.6 1.1 ± 0.1 6.8 ± 0.7 0.32 455/337

χ = 60◦ , ξ = 10◦
Observation ID a
N (10 )23 b
K (10 ) −6
α b
βb Jc (10−3 ) T in d χ2 /dof
(keV/cm2 /s) (keV/cm2 /s) (keV) (keV/cm2 /s) (keV)

0693760401 4.8 ± 0.3 0.50 ± 0.35 -7.7 ± 0.9 0.83 ± 0.09 4.7 ± 0.3 0.32 334/197
0748390901 2.0 ± 0.1 1.6 ± 0.9 -6.6 ± 0.6 1.1 ± 0.1 7.1 ± 0.7 0.32 457/337

Notes.
(a)
torusdisk norm
(b)
cutoffpl: norm, photon index, and e-folding energy of exponential cut-off
(c)
bbody norm
(d)
Temperature at the inner radius of the accretion disk

and upper limit, χ = 20◦ , ξ = 10◦ and χ = 60◦ , ξ = 10◦ , respec- be adjusted in the fit. In addition, Rodríguez Castillo et al. (2020)
tively. This is not surprising since the source is very variable and also adopted the absorption model tbabs, an updated version of
its PF does not increase monotonically with the energy, as shown wabs, and this can explain the higher NH .
in Figure 3 in Brightman et al. (2022). The reasons behind the For a magnetic field B = 8 × 1012 G, we fitted the data
source variability are not well known. They could be related to with the superposition of three components: a phenomenological
the fact that the source entered a propeller stage, or because of a blackbody, the torusdisk model, and a phenomenological cutoff
variation in the mass transfer rate, which is unstable. power law, which represents the non-thermal emission from the
We used both the geometries to perform the spectral com- accretion column. With the three component model, the spec-
parison, adding to the torusdiks multicolor blackbody the phe- trum is reproduced with great accuracy (Figure 7). A similar 3-
nomenological components to fit the entire spectrum. We consid- component spectral analysis was also performed by Brightman
ered two representative values for the magnetic field, B = 1012 et al. (2022). They fitted the data using the same blackbody and
G and 8 × 1012 G, respectively, which brackets the true value. cutoff power law components, and a phenomenological multi-
For each of them, we computed the Alfvén radius rA and the color blackbody at intermediate energies, well reproducing the
torus temperature. We also fixed the disk outer radius, Rdisk = whole spectrum. For the same reasons repeated above we did
200 RNS , and varied the temperature at the inner radius of the not use their best-fit parameters for the blackbody and the cut-
disk between 0.25 and 0.5 keV. off power law, but we left them free to vary (apart from those
that we fixed to ensure the convergence of the fit: the bbody
For a magnetic field B = 1012 G we were able to fit the
temperature and the total column density). The T in is ∼ 0.36
spectrum with two thermal components: a blackbody at lower
keV, compatible with our value, while the column density NH
energies and our torusdik model at higher energies (figure 6),
is lower than that found by Rodríguez Castillo et al. (2020),
similar to what was done by Rodríguez Castillo et al. (2020).
NH ∼ 3.3 × 1020 ± 0.7 cm−2 , and closer to our value (the residual
In the fitting procedure, we left free to vary the normalization
difference being again possibly related to the different absorp-
of the torusdisk model and the column density NH , while we
tion model). For the cutoff power law parameters, Brightman
fixed at 0.25 keV the bbody temperature, and we let the tem-
et al. (2022) fixed them using the values found by Walton et al.
perature T in to vary between 0.25 and 0.5 keV, but fixing it for
(2018a), Ecut (here β) ∼ 8.1 keV and Γ(here α) ∼ 0.8, that are
each simulated spectrum. We did not fix the parameters of the
different from our results.
phenomenological bbody component to the values inferred by
Rodríguez Castillo et al. (2020), because our multicolor black- From spectral comparison we obtained good constraints on
body is different. The two thermal component models of Ro- the temperature T in and on the magnetic field strength B, but no
dríguez Castillo et al. (2020) fit well the spectrum, with a higher additional constraints on the geometry of view. The latter can be
column density of 5.9 × 1020 cm−2 , and a higher value for the further investigated through polarization analyses.
temperature and radius of the hotter thermal component (1.5 keV
and 97 RNS , respectively), while T in is compatible within the un- 5.1.2. Polarization
certainties with our value. On the other hand, our physically mo-
tivated model (blackbody + torusdisk) does not provide a sat- As shown in Sect. 4.3.1 the polarization degree increases with
isfactory fit to the data (figure 6) because of residuals between the energy (figure 9). This is a consequence of the strong mag-
∼ 1.5 and 3.0 keV. The reason is that in our torusdisk multicolor netic field considered. For B = 8 × 1012 G, high energy photons
blackbody, the energy range of the torus temperature varies de- are emitted from the torus inside the adiabatic region. This ex-
pending on the choice of the magnetic field strength and cannot plains the oscillating behavior of the polarization degree with
Article number, page 13 of 17
 A&A proofs: manuscript no. aanda

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

the phase, which is simply a consequence of the rotation of larization degree remains the same even changing the intrinsic
the source combined with the geometry of view, that changes polarization fraction ΠL .
the emitting regions in view from which the polarized radiation The polarization angle shows an oscillating behavior with the
comes. For regions closer to the equator PD is low because pho- phase, going from 90◦ to 0◦ or 180◦ , as expected from a rotating
tons have lower energies (the torus temperature is lower) and are vector model, while remaining almost constant with the energy.
emitted outside the adiabatic region. The oscillations are caused by the geometry of view combined
with the rotation of the star. Because the part in view of the emit-
Fixing the energy and the phase, the polarization degree in- ting components varies from regions emitting polarized radiation
creases with the angle χ (fig. 9). This is because of the geomet- (X-mode photons), where we expect the polarization angle to os-
rical depolarization effect. To understand this point, let us con- cillate around 90◦ , to regions emitting mainly unpolarized radia-
sider two limiting cases: a pole-on configuration (χ = ξ) and an tion, that give a PA of 0◦ or 180◦ , if we recall how PA is defined
equator-on configuration (χ − ξ = 90◦ ), with respect to the mag- (equation (9)), and that the values of U and Q, in those regions,
netic axis. In the first case, there is no preferential direction of are close to zero.
the magnetic field vector, since the field assumes a radial config-
uration. Therefore photons escaping from those regions (linearly
polarized parallel to or perpendicular to the plane of the vectors 5.2. NGC 7793 P13
k and B) will have an electric field with a radial configuration, 5.2.1. Light curves and spectra
which means to have opposite polarization contribution, so a null
total polarization. In the second case, there is a preferential di- We repeated the previous analysis on NGC 7793 P13, start-
rection of the magnetic field vector, so summing the polarization ing from the light curve comparisons. We found two limit-
contributions we have a non-zero polarization degree. For the ing viewing geometries (this time with a smaller range of χ):
χ = 20◦ , ξ = 10◦ geometry, we are closer to a pole-on con- χ = 40◦ , ξ = 10◦ , and χ = 60◦ , ξ = 10◦ . The pulsed frac-
figuration, while increasing χ up to 60◦ we move towards an tion of NGC 7793 P13 presents a more regular behavior, in-
equator-on configuration. creasing quasi-monotonically with energy (Israel et al. 2017b).
We used both geometries to perform the spectral fits, consid-
It is also worth noting that even if we cannot distinguish the ering a magnetic field strength of 4 × 1012 G (computing the
viewing geometry from the X-ray spectrum, polarization mea- torus temperature, the Alfvén radius, table 3, and range of T in
surements may allow us to do it. Indeed, the behavior of the po- as before), and setting the disk outer radius Rdisk = 200 RNS .
Article number, page 14 of 17
 S. Conforti et al.: Pulsating ultraluminous X-ray sources: modeling the thermal emission and polarization properties

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

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◦ .

The spectrum is well reproduced by the superposition of a X-ray energy band is thermal and well described by a multicolor
blackbody, the torusdisk model, and a cutoff power law. Israel blackbody.
et al. (2017b) performed a phenomenological fit using an ab-
sorbed power law with a high energy cut-off and a multicolor
blackbody (phabs[highecut*(powerlaw+bbodyrad)]), obtaining 5.2.2. Polarization
the following parameters: NH = 9.6×1020 cm−2 , Γ(here α) ∼ 1.2,
Ecut (here β) ∼ 6 keV, kT BB ∼ 0.2 keV. Despite some differences As for M51 ULX-7, also in NGC 7793 P13 for a given geometry
in the fitting parameters, both models agree on the type of emis- of view the polarization degree increases with energy, because
sion at high energies, well described by a power law, probably of the strong magnetic field, and it oscillates with phase. The
produced by an accretion column, while the emission in the soft maximum polarization degree is around 4.7%, lower than that of
M51 ULX-7 (∼ 19%), using the lower magnetic field strength,
Article number, page 15 of 17
 A&A proofs: manuscript no. aanda

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◦ .

B = 4 × 1012 G. For such magnetic field we expect the radiation and ∼ 1013 G. The accreting material originates from the donor
to be poorly polarized (in the X-mode), expecting less variation star and proceeds through a geometrically-thin accretion disk up
of PD between more energetic and less energetic radiation (fig- to the Alfvén radius, where particles are funneled along the mag-
ure 10, top right). We observe oscillations with the phase still re- netic field lines, which follow a dipolar field topology, towards
lated to the rotation of the source combined with the geometry of the magnetic poles of the star. Each point of the torus and the
view. The reason why this time the polarization degree decreases disk emits like a blackbody, with a local temperature that de-
with χ depends on the different magnetic field strength used, on pends on the magnetic colatitude θ, if the point belongs to the
the geometry of view, and on how the geometrical depolariza- torus, or on the radial distance if the point belongs to the disk.
tion effect dominates over the vacuum birefringence and vice The emission of the whole source is then a complex multicolor
versa. Indeed, despite the favorable viewing angle (χ = 60◦ ), blackbody. We considered as well polarization observables, tak-
because of the depolarization effect mentioned above, the lower ing as polarized only the radiation coming from below the adia-
magnetic field strength reduces the polarized radiation. As a re- batic radius rpl , while that coming from the disk and the part of
sult, in regions near the equator the radiation is so weakly po- the torus outside the region of adiabatic propagation is unpolar-
larized that, even summing all the polarization contributions, we ized. We computed the flux and the polarization observables for
cannot achieve a higher polarization degree (figure 10, top left). different viewing geometries (χ, ξ) of each source.
Also for NGC 7793 P13, the behavior of PD may allow us to We made direct comparisons of the light curves considering
distinguish the geometry of view if observed polarization mea- the 2–3 keV and 3–4 keV energy bands. We compared each ob-
surements are available. served profile with six simulations having different geometries
The polarization angle oscillates with phase and remains of view, and from these comparisons, we constrained the view-
constant with energy, but switching drastically between the high- ing geometry obtaining a lower and upper limit of the angle χ,
est (∼ 180◦ ) and lowest value (∼ 1◦ ) at lower energies for the which was then used in the spectral analysis.
geometry of view χ = 60◦ ξ = 10◦ . The reason is probably that We compared spectra from different observations with those
the radiation is mainly unpolarized. 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
6. Conclusions observations (0824450901, 0830191501 and 0830191601), and
In this work, we presented a simplified model that reproduces compared the spectrum with a multicomponent model that com-
the X-ray thermal emission of pulsating ultraluminous X-ray bines a blackbody (bbody, between 0.1 − 1.5 keV), a multicolor
sources, in terms of emission from an accretion disk and an ac- blackbody (output of the model, between 1.5 − 5 keV), and a
cretion envelope modeled as a torus confined by the magnetic cut-off power law (cutoffpl, between 5 − 10 keV), assuming a
field lines and extending up to the magnetospheric radius. We magnetic field strength B = 8 × 1012 G. The best fitting internal
used and modified the ray-tracer code discussed in Taverna et al. disk temperature is T in ≈ 0.3 keV. We did the same for NGC
(2015), to reproduce spectra and light curves of PULXs in two 7793 P13, using the XMM-Newton observations 0693760401
different energy bands. The model also simulates the polariza- and 0748390901, finding a T in ≈ 0.32 keV.
tion degree and angle of the radiation emitted as a function of For a given viewing geometry, the polarization degree in-
energy and phase. We tested the model on two PULXs: M51 creases with the energy and shows an oscillating behavior with
ULX-7 and NGC 7793 P13, using for both XMM-Newton obser- the phase. The former is a consequence of the strong magnetic
vations. The aim was to test the validity of our model and derive field considered in the model, while the oscillating behavior with
information on the geometry of view, magnetic field, thermal and the phase is a consequence of the rotation of the source combined
polarization properties of such sources. with the geometry of view. Fixing the energy and phase but vary-
In the model we considered the emission from a highly mag- ing the angle χ we note that the value of the polarization degree
netized neutron star in an accreting binary system with a mass of increases. The variation of PD from one geometry of view to an-
1.4 M⊙ , radius 10 km, and magnetic field strength between 1012 other could allow us to distinguish them if observed polarization
Article number, page 16 of 17
 S. Conforti et al.: Pulsating ultraluminous X-ray sources: modeling the thermal emission and polarization properties

measurements are available. Therefore, the polarization analysis Middleton, M. J., Walton, D. J., Fabian, A., et al. 2015b, Monthly Notices of the
has proved to be crucial in determining the geometry of view of Royal Astronomical Society, 454, 3134
the source when the spectral one gives degenerate results. Miller, J., Fabbiano, G., Miller, M., & Fabian, A. 2003, The Astrophysical Jour-
nal, 585, L37
The model reproduces well the observed spectra of the two Morrison, R. & McCammon, D. 1983, Astrophysical Journal, Part 1 (ISSN 0004-
pulsating ULXs, particularly in the 1.5–5 keV interval, where 637X), vol. 270, July 1, 1983, p. 119-122., 270, 119
the contribution of the torus with the disk is dominant. Motch, C., Pakull, M., Grisé, F., & Soria, R. 2011, Astronomische Nachrichten,
332, 367
Acknowledgements. SC is supported by the INAF Doctorate program. SC and Mushtukov, A. A., Suleimanov, V. F., Tsygankov, S. S., & Ingram, A. 2017,
LZ acknowledge financial support from the INAF Research Grant “Uncover- Monthly Notices of the Royal Astronomical Society, 467, 1202
ing the optical beat of the fastest magnetized neutron stars (FANS)". RTa and Mushtukov, A. A. et al. 2015, Monthly Notices of the Royal Astronomical Soci-
RTu acknowledge financial support from the Italian MUR through grant PRIN ety, 454, 2539
2022LWPEXW. SC thanks E. Giunchi, B. Ambrosio and I. Salmaso for technical Pinto, C., Mehdipour, M., Walton, D. J., et al. 2020, Monthly Notices of the
support. Royal Astronomical Society, 491, 5702
Pinto, C., Middleton, M. J., & Fabian, A. C. 2016, Nature, 533, 64
Pintore, F. & Zampieri, L. 2012, Monthly Notices of the Royal Astronomical
Society, 420, 1107
References Poutanen, J., Lipunova, G., Fabrika, S., Butkevich, A. G., & Abolmasov, P.
2007a, Monthly Notices of the Royal Astronomical Society, 377, 1187
Abolmasov, P., Fabrika, S., Sholukhova, O., & Afanasiev, V. 2007, Astrophysical Poutanen, J. et al. 2007b, Monthly Notices of the Royal Astronomical Society,
Bulletin, 62, 36 377, 1187
Arnaud, K. 1996, in ASP Conf., Vol. 17 Quintin, E., Webb, N. A., Gúrpide, A., Bachetti, M., & Fürst, F. 2021, Monthly
Bachetti et al. 2014, Nature, 514, 202 Notices of the Royal Astronomical Society, 503, 5485
Basko, M. & Sunyaev, R. A. 1976, Monthly Notices of the Royal Astronomical Roberts, T. & Warwick, R. 2000, Monthly Notices of the Royal Astronomical
Society, 175, 395 Society, 315, 98
Brice, N., Zane, S., Taverna, R., Turolla, R., & Wu, K. 2023, Monthly Notices of Rodríguez Castillo, G., Israel, G. L., Belfiore, A., et al. 2020, The Astrophysical
the Royal Astronomical Society, 525, 4176 Journal, 895, 60
Brice, N., Zane, S., Turolla, R., & Wu, K. 2021, Monthly Notices of the Royal Rybicki, G. B. & Lightman, A. P. 1991, Radiative processes in astrophysics (John
Astronomical Society, 504, 701 Wiley & Sons)
Brightman, M., Bachetti, M., Earnshaw, H., et al. 2022, The Astrophysical Jour- Sathyaprakash, R. et al. 2019, Monthly Notices of the Royal Astronomical Soci-
nal, 925, 18 ety: Letters, 488, L35
Carpano, S. et al. 2018, Monthly Notices of the Royal Astronomical Society: Shakura, N. I. & Sunyaev, R. A. 1973, Astronomy and Astrophysics, Vol. 24, p.
Letters, 476, L45 337-355, 24, 337
Colbert, E. J. & Mushotzky, R. F. 1999, The Astrophysical Journal, 519, 89 Stobbart, A.-M., Roberts, T., & Wilms, J. 2006, Monthly Notices of the Royal
Fabbiano, G. 1989, Annual review of astronomy and astrophysics, 27, 87 Astronomical Society, 368, 397
Fabbiano, G., Kim, D.-W., & Trinchieri, G. 1992, The Astrophysical Journal Sutton, A. D., Roberts, T. P., Gladstone, J. C., et al. 2013a, Monthly Notices of
Supplement Series, 80, 531 the Royal Astronomical Society, 434, 1702
Fabrika, S. & Mescheryakov, A. 2001, in Symposium-International Astronomi- Sutton, A. D., Roberts, T. P., & Middleton, M. J. 2013b, Monthly Notices of the
cal Union, Vol. 205, Cambridge University Press, 268–269 Royal Astronomical Society, 435, 1758
Feng, H. & Soria, R. 2011, New Astronomy Reviews, 55, 166 Swartz, D. A. et al. 2004, The Astrophysical Journal Supplement Series, 154,
Fernández, R. & Davis, S. W. 2011, The Astrophysical Journal, 730, 131 519
Fürst, F., Walton, D., Harrison, F., et al. 2016, The Astrophysical Journal Letters, Taverna, R., Muleri, F., Turolla, R., et al. 2014, Monthly Notices of the Royal
831, L14 Astronomical Society, 438, 1686
Fürst, F., Walton, D., Heida, M., et al. 2021, Astronomy & Astrophysics, 651, Taverna, R. & Turolla, R. 2017, Monthly Notices of the Royal Astronomical
A75 Society, 469, 3610
Gladstone, J. C., Roberts, T. P., & Done, C. 2009, Monthly Notices of the Royal Taverna, R., Turolla, R., Gonzalez Caniulef, D., et al. 2015, Monthly Notices of
Astronomical Society, 397, 1836 the Royal Astronomical Society, 454, 3254
Gnedin, Y. N. & Pavlov, G. 1974, Soviet Journal of Experimental and Theoretical Terashima, Y. & Wilson, A. S. 2004, The Astrophysical Journal, 601, 735
Physics, 38, 903 Vasilopoulos, G., Lander, S., Koliopanos, F., & Bailyn, C. 2020, Monthly Notices
Harding, A. K. & Lai, D. 2006, Reports on Progress in Physics, 69, 2631 of the Royal Astronomical Society, 491, 4949
Heisenberg, W. & Euler, H. 1936, Zeitschrift für Physik, 98, 714 Ventura, J. 1979, Physical Review D, 19, 1684
Herold, H. 1979, Physical Review D, 19, 2868 Walton, D. 2015, XMM-Newton Proposal, 68
Heyl, J. S. & Shaviv, N. J. 2002, Physical Review D, 66, 023002 Walton, D., Fürst, F., Harrison, F., et al. 2018a, Monthly Notices of the Royal
Heyl, J. S., Shaviv, N. J., & Lloyd, D. 2003, Monthly Notices of the Royal As- Astronomical Society, 473, 4360
tronomical Society, 342, 134 Walton, D., Fürst, F., Heida, M., et al. 2018b, The Astrophysical Journal, 856,
Ho, W. C. & Lai, D. 2003, Monthly Notices of the Royal Astronomical Society, 128
338, 233 Weisskopf, M. C., Soffitta, P., Baldini, L., et al. 2022, Journal of Astronomical
Israel, G. L., Belfiore, A., Stella, L., et al. 2017a, Science, 355, 817 Telescopes, Instruments, and Systems, 8, 026002
Israel, G. L., Papitto, A., Esposito, P., et al. 2017b, Monthly Notices of the Royal Zampieri, L. & Roberts, T. 2009, Monthly Notices of The Royal Astronomical
Astronomical Society: Letters, 466, L48 Society, 400, 677
Kaaret, P., Feng, H., & Roberts, T. P. 2017, Annual Review of Astronomy and
Astrophysics, 55, 303
Kaaret, P. et al. 2001, Monthly Notices of the Royal Astronomical Society, 321,
L29
King, A. R. et al. 2001, The Astrophysical Journal, 552, L109
Lai, D., Ho, W., van Adelsberg, M., Wang, C., & Heyl, J. 2010, X-ray Polarime-
try: A New Window in Astrophysics, 157
Liu, Q. & Mirabel, I. 2005, Astronomy & Astrophysics, 429, 1125
Long, K. S., Dodorico, S., Charles, P. A., & Dopita, M. A. 1981, The Astrophys-
ical Journal, 246, L61
Lyubarskii, Y. E. & Syunyaev, R. 1988, Soviet Astronomy Letters, 14, 390
Makishima, K. et al. 2000, The Astrophysical Journal, 535, 632
Mészáros, P. 1992, High-energy radiation from magnetized neutron stars (Uni-
versity of Chicago press)
Middleton, M. J., Heil, L., Pintore, F., Walton, D. J., & Roberts, T. P. 2015a,
Monthly Notices of the Royal Astronomical Society, 447, 3243
Middleton, M. J., Higginbottom, N., Knigge, C., Khan, N., & Wiktorowicz, G.
2022, Monthly Notices of the Royal Astronomical Society, 509, 1119

Article number, page 17 of 17