© The Authors 2025

1. Introduction

Ultraluminous X-ray sources (ULXs) provide valuable insights into super-Eddington accretion onto compact objects located in nearby galaxies (Kaaret et al. 2017). The nature of these point sources, with bolometric luminosities equal to or greater than 10³⁹ erg s⁻¹, remains unclear. Several hypotheses attempt to explain the observed high X-ray luminosities. One suggests that this phenomenon can be produced by the accretion of intermediate-mass black holes (IMBHs) at sub-Eddington rates (Li & van den Heuvel 1999; Madhusudhan et al. 2006; Strohmayer & Mushotzky 2009). Another proposes that compact stellar-mass objects with jets undergo sub-Eddington accretion regimes (Begelman et al. 2006; Saavedra et al. 2023a). Additionally, stellar-mass black holes (BHs) or neutron stars (NSs) with supercritical accretion rates have been proposed (Poutanen et al. 2007; Koliopanos et al. 2017; Mushtukov et al. 2017; Abarca et al. 2018; Fabrika et al. 2021; Mills et al. 2024).

Ultraluminous X-ray sources are classified into distinct states based on their spectral components and hardness, as documented in previous studies (Sutton et al. 2013; Kaaret et al. 2017; Fabrika et al. 2021). The dominant state is the ultraluminous state, characterized by a cool accretion disk (kT < 0.5 keV) and a power-law tail. Within this state, ULXs can be further classified as hard ultraluminous (HUL) or soft ultraluminous (SUL) sources, depending on their spectral hardness (Sutton et al. 2013). The hard intermediate and soft bright states are related classifications of these two regimes. In soft bright states, ULXs typically exhibit higher luminosities compared to hard intermediate states, due to the down-scattering of photons by the dense medium of clumpy winds at elevated accretion rates. Other ULX states include the broadened disk (BD) state and the super-soft ultraluminous (SSUL) state. BD states are characterized by thermal emission from what seems to be a geometrically modified accretion disk, while SSUL states are dominated by a single-component cool blackbody component. It is important to note that the boundary between these classifications is not always clear, with a significant fraction of ULXs exhibiting spectral transitions between different states (e.g., Walton et al. 2020; D’Aì et al. 2021). These transitions are likely the result of variations in accretion rates, disk occultations, wind outflow strengths, and several other physical parameters.

In recent years, observations with advanced X-ray telescopes such as XMM-Newton and the Nuclear Spectroscopic Telescope Array (NuSTAR) have facilitated the identification of the majority of ULX sources as super-Eddington accreting X-ray binaries (Walton et al. 2013; Mukherjee et al. 2015). This finding was further strengthened by the remarkable discovery that several ULXs are actually powered by accreting NSs (Bachetti et al. 2014), the so-called pulsating ULXs (PULXs).

The first conceptualizations of super-Eddington accretion involved a geometrically thin and optically thick disk surrounding the central BH (Shakura & Sunyaev 1973). A key prediction of these models is the generation of substantial radiation-driven outflows of ionized gas, emanating from regions proximal to the BH and reaching slightly relativistic velocities. However, current models mainly consider disk inflation within the critical radius, where the radiation pressure exceeds the gravitational attraction. Within this inner region, the disk becomes unstable, and a dense, opaque wind evacuates most of the matter. Photon-trapping and advection extract energy from the disk and feed the BH at a rate close to the Eddington limit (Ohsuga et al. 2003; Fukue 2004; Dotan & Shaviv 2011; Ogawa et al. 2017).

ESO 501–023 is a barred spiral galaxy (SBdm; de Vaucouleurs et al. 1991) at a distance of ∼7 Mpc (Karachentsev et al. 2013). Recent Chandra observations have revealed a bright ULX present in ESO 501–023 (Somers et al. 2013). The spectral properties of ESO 501–023 ULX (hereafter ESO501) show a photon index of Γ ∼ 1.2 and an unabsorbed flux of ∼1.1 × 10⁻¹² erg s⁻¹ cm⁻², which corresponds to a luminosity of ∼5 × 10³⁹ erg s⁻¹ for a distance of ∼7 Mpc (Somers et al. 2013) in the 0.3–8 keV energy range. We note, however, that Somers et al. (2013) suggested that ESO501 is located farther away on the basis of spectral/luminosity considerations (Sutton et al. 2013), and so the luminosity remains unclear.

IC 5052 is a barred spiral galaxy (SBd; de Vaucouleurs et al. 1991) located closer than ESO 501–023 at a distance of about ∼5.5 Mpc (Tully et al. 2013). It is classified as a late-type galaxy with a low star formation rate (taken from the NED). During a survey of X-ray sources in nearby late-type galaxies, Chatterjee et al. (2015) discovered a ULX in IC 5052. This ULX, designated IC 5052 X-1 (hereafter IC5052), has a luminosity of ∼2.9 × 10³⁹ erg s⁻¹. The spectral properties of IC 5052 X-1 include a photon index of ∼1.7 and a disk blackbody temperature of ∼0.45 keV in the 0.2–10 keV energy range (Chatterjee et al. 2015). They suggest that IC5052 may be an IMBH based on the low accretion disk temperatures, which fit the bimodal temperature distribution found by Winter et al. (2006).

In this study we analyzed archival XMM-Newton and NuSTAR observations of the ULX sources in the galaxies ESO 501–023 and IC 5052. Our goal is to refine the understanding of the geometry and structure of the accretion disk around the compact objects, to identify spectral hints of outflows, and to explore possible pulsations.

The paper is organized as follows. In Sect. 2 we describe the observations and the data reduction methods used in this study, detailing the processing of the XMM-Newton and NuSTAR data. In Sect. 3 we present the main results of our timing and spectral analyses, including the investigation of possible pulsations and the characterization of the spectral states of ESO501 and IC5052. In Sect. 4 we provide a detailed discussion of the physical interpretation of these results, focusing on the modeling of super-Eddington accretion and the possible geometries of the accretion disks. Finally, Sect. 5 summarizes our main findings and discusses their broader implications for the study of ULX sources.

2. Observations and data analysis

2.1. XMM-Newton data

XMM-Newton performed observations of IC5052 on 5 October 2022 with an exposure time of ∼61.9 ks (ObsID 0912390101) and of ESO501 on 20 November 2022 with an exposure time of ∼53 ks (ObsID 0912390201). During both observations, the MOS and PN cameras were set to full-frame mode.

XMM-Newton data were processed using version 21.0.0 of the Science Analysis System (SAS) software together with calibration and configuration files available until October 2023.We generated calibrated event files by filtering out high-energy background activity using standard count rate limits for PN and MOS instruments. We obtained good time intervals of ∼21 ks (PN camera) for IC5052 and ∼39 ks (PN camera) for ESO501. Source spectra and light curves were extracted using circular regions with a radius of 25 arcsec and 35 arcsec for IC5052 and ESO501, respectively. Background regions were positioned on the same chip at a comparable distance from the readout node. Event selection criteria included FLAG==0, with PATTERN ≤ 12 for the EPIC-MOS cameras and PATTERN ≤ 4 for the EPIC-PN cameras. In all cases the spectra were binned to ensure a minimum of 20 counts per bin, with a oversample==3 in the 0.3–12 keV energy range to allow χ² statistics to be applied. Redistribution matrices and auxiliary response files were generated using the rmfgen and arfgen tasks, respectively.

Light curves in the 0.3–12 keV energy range were corrected for background and exposure using the epiclccorr task of SAS. Consequently, a binning of 140 seconds was used to visualize the light curves. Barycentric corrections were applied using the SAS barycen task.

2.2. NuSTAR data

NuSTAR (Harrison et al. 2013) is an X-ray satellite equipped with two detectors, FPMA and FPMB, operating in the 3 to 79 keV energy band. NuSTAR observed IC5052 in October 2022 (ObsID 30801029002) with an exposure of ∼113 ks and ESO501 in November 2022 (ObsID 30801030002) with an exposure of ∼107 ks. We analyzed the data using the NuSTARDAS-v. 2.0.0 analysis software from HEASoft v.6.32.1 and calibration files from CALDB (v.20230816). We took source events that accumulated within a circular region of 60 arcsec radius centered on the point spread function. Background regions for the FPMA/B modules were selected far from source and detector edge contamination with a radius of 60 arcsec.

We used the nupipeline task to create level 2 calibrated event files, with saacalc=2, saamode=OPTIMIZED, tentacle=YES to filter for South Atlantic Anomaly passing epochs for ESO501 and with saacalc=1, saamode=OPTIMIZED, tentacle=YES for IC5052, both obtained from the South Atlantic Anomaly filtering report¹. This resulted in a loss of ∼0.8% of the total unfiltered exposure time. We extracted light curves and spectra using the nuproducts task. We applied barycentric correction using the barycorr task with the nuCclock20100101v173 clock correction file. We used the celestial coordinates α = 158.8425 deg and δ = −24.75358 deg for ESO501 and α = 313.0706 and δ = −69.22131 for IC5052, with JPL-DE200 solar system ephemeris. We produced background corrected light curves for each module, and subsequently merged the two light curves using the LCMATH task. The spectra were binned to ensure a minimum of 20 counts per energy bin in the 3–79 keV energy range, allowing the χ² statistical test to be applied.

3. Results

3.1. Timing analysis

We analyzed the background-subtracted light curves for both ULXs using combined data from the XMM-Newton and NuSTAR observations to search for any significant variability. The light curves were extracted in the 0.3–12 keV energy range for XMM-Newton and in the 3–78 keV range for NuSTAR. For each instrument, count rates were calculated to assess the presence of any high or low energy variability across the observing time.

For the XMM-Newton observations, we found no significant variability over the entire observation time for either source. For the PN camera, the count rate for ESO501 was 0.432 ± 0.003 cts s⁻¹ (with an exposure time of 39.2 ks), for MOS1 it was 0.117 ± 0.002 cts s⁻¹ (48 ks), and for MOS2 it was 0.123 ± 0.002 cts s⁻¹ (47.8 ks). For IC5052, the count rate for the PN camera was 0.162 ± 0.003 cts s⁻¹ (21 ks), for MOS1 it was 0.039 ± 0.001 cts s⁻¹ (47.6 ks), and for MOS2 it was 0.053 ± 0.001 cts s⁻¹ (48.4 ks).

Similarly, the NuSTAR observations showed no significant variability for either source. For FPMA on ESO501, the count rate was 0.0096 ± 0.0004 cts s⁻¹ (102.1 ks) and for FPMB, the count rate was 0.0084 ± 0.0005 cts s⁻¹ (101.3 ks). For IC5052, the count rate for FPMA was 0.0099 ± 0.0004 cts s⁻¹ (106.6 ks) and for FPMB it was 0.0093 ± 0.0004 cts s⁻¹ (105.2 ks).

We used the spectral timing algorithms provided by the HENDRICS (Bachetti 2018) and Stingray (Huppenkothen et al. 2019) software packages to investigate potential pulsations within the XMM-Newton and NuSTAR calibrated event files, using the same energy ranges as above. We performed this search in frequency ranges appropriate for each instrument, based on their sampling times and observation mode.

For the XMM-Newton EPIC-PN data, we searched for pulsations in the 0.01–6 Hz frequency range since the full-frame mode provides a Nyquist frequency of about 6 Hz. For the EPIC-MOS1/2 cameras, due to their longer readout times, we restricted the frequency search to the 0.01–0.2 Hz range, which corresponds to their Nyquist frequencies in full-frame mode. For NuSTAR, despite the shorter readout times of the FPMA/B detectors, which allow a theoretical Nyquist frequency of up 200 Hz, we restricted our search to the 0.01–10.0 Hz range to avoid high frequency noise and to be consistent with the XMM-Newton analysis.

We used the HENaccelsearch tool to search for pulsation candidates in the 0.01–10 Hz range, allowing for variable period derivatives up to 10⁻⁹ Hz s⁻¹. No significant pulsation candidates were detected in any of the datasets. We then applied the HENz2vspf tool to determine an upper bound on the pulse fraction (PF) potentially exhibited by each signal. This procedure scrambles the event times and adds a pulsation with a random pulsed fraction. It then extracts the maximum of the Z² distribution within a small interval around the pulsation. This process is repeated N_trial times. The PF is defined as the ratio between the amplitude of the pulsation and the maximum flux. Based on 1000 simulations, we derived an upper limit for the PF (with 90% confidence). Specifically, for ESO501, the PF upper limit was ∼9% for the XMM-Newton data and ∼7% for the NuSTAR data. Similarly, for IC5052, the upper limit was ∼11% for XMM-Newton and ∼6% for NuSTAR.

3.2. Spectral analysis

The spectral analysis of the X-ray observations was performed using XSPEC v.12.13.1 (Arnaud 1996). Unabsorbed fluxes and luminosities were estimated using the cflux and cglumin convolution models, assuming distances of 7 Mpc and 5.5 Mpc for ESO 501-023 and IC 5052, respectively. We used Markov chain Monte Carlo (MCMC) chains to estimate uncertainties by generating 10⁶ samples using XSPEC’s chain command. The number of walkers was set to 144, corresponding to 12 times the number of free parameters (12 in total). This choice was motivated by two factors: first, having a sufficient number of walkers to ensure a thorough exploration of the parameter space, as having multiple walkers helps avoid chains becoming trapped in local minima; and second, the availability of 12 processors allowed us to parallelize the computation across all cores. This setup optimized the performance of the MCMC process and significantly accelerated convergence through parallel processing. We ensured successful chain convergence by inspecting each parameter trace plot and verifying that sufficient state changes were reached by all walkers (i.e., autocorrelation time close to unity; see Fogantini et al. 2023 and Saavedra et al. 2023b for more details).

In order to fully characterize the broadband average spectra of the two sources, the five instruments (PN + MOS1/2 and FPMA/B) were used simultaneously. For ESO501, the source was significantly detected above background between 0.3 and 10 eV in the XMM-Newton data and between 3 and 20 keV in the NuSTAR data. On the same note, for IC5052, the source is above background between 0.5 and 20 keV in the XMM-Newton + NuSTAR spectra. We included a calibration constant (constant model in XSPEC) to account for systematic differences between instruments.

Following Sutton et al. (2013), we began by fitting the average spectra of the two sources in the 0.3–10.0 keV band with a model that includes a multicolor disk (diskbb) and a power-law component. For ESO501 we obtained a disk temperature of kT_in = 1.37 ± 0.06 keV and a photon index of Γ = 2.7 ± 0.3, with a reduced χ² value of ∼1.03 (407 d.o.f.). The unabsorbed flux ratio of F_PL/F_disk in the 0.3–10 keV band was $4.7_{-0.9}^{+1.6}$. For IC 5052 we obtained a disk temperature of kT${}_{\mathrm{in}}=0.{21}_{-0.03}^{+0.05}$ keV and a photon index of Γ = 1.82 ± 0.10, with a reduced χ² value of ∼0.98 (332 d.o.f.). Using the decision tree described in Fig. 2 of Sutton et al. (2013), and based on comparing our results in Fig. 1 with those in Fig. 1 of Sutton et al. (2013), we can classify ESO501 and IC5052 into the BD and HUL states, respectively. Additionally, the emission peak observed around 5 keV in the spectrum of ESO501 is consistent with the BD state (see Fig. 1).

Fig. 1. XMM-Newton+NuSTAR unfolded spectra and residuals of ESO501 (left column) and IC5052 (right column). Each model and best-fit statistic are specified.

We formulated two preliminary hypotheses to investigate the properties of ESO501 and IC5052. The first hypothesis was for a nonmagnetic accretor, such as a BH, while the second hypothesis was for a magnetic accretor, such as a PULX, despite the absence of detected pulsations in both sources. This study included two neutral absorption models, the first fitted to galactic absorption columns of 4.6 × 10²⁰ cm⁻² and 4 × 10²⁰ cm⁻² for ESO501 and IC5052, respectively, derived from the NH web tool². The second TBABS model was left free to account for local absorption. For the neutral absorption we used abundances from Wilms et al. (2000) and cross sections from Verner et al. (1996).

For nonmagnetic objects, a two-component thermal model was used to phenomenologically describe the accretion flow scenario. The DISKBB model (Mitsuda et al. 1984) was used to characterize the accretion flow in the outer region, whereas the DISKPBB model (Mineshige et al. 1994) was used to describe the flow in the inner regions where a departure from the thin disk approximation is expected. While the DISKPBB model assumes a thin disk profile with a fixed exponent p value of 0.75, the DISKPBB model allows p to vary as a free parameter. This model combination is often used to fit low energy X-ray data (E < 10 keV) when studying ULX spectra. The emission at high energies (E > 10 keV) is thought to result from Compton scattering of disk photons in a hot electron plasma. To avoid inaccurate extrapolations of this emission at low energies, the SIMPL model (Steiner et al. 2009) is used. This complementary component works with the diskpbb model, which characterizes the disk’s hottest region. The full model in XSPEC reads tbabs×tbabs×(diskbb+simpl×diskpbb), which we refer to this model as simpl hereafter.

The description of a magnetic accretor, specifically an PULX, was derived from the framework presented in Walton et al. (2018c,b). This framework includes two thermal blackbody components, representing the accretion flow from the magnetosphere, and a power-law component with an exponential cutoff (cutoffpl), which delineates the central accretion columns produced by the material flow at the magnetic poles. The thermal component is modeled using the diskbb+diskpbb combination commonly used in PULX studies. Parameters for the high-energy part of the spectrum (E > 10 keV) were determined from averages derived from the pulsed emission of known PULX. These parameters included Γ = 0.59 and E_cut = 7.9 keV for the cutoffpl component (Brightman et al. 2016; Walton et al. 2018c). The full model in XSPEC reads tbabs×tbabs×(diskbb+diskpbb+cutoffpl), which we call cutoffpl hereafter for simplicity. The results of applying these models to the XMM-Newton+NuSTAR data of both sources are shown in Table 1 and Fig. 1.

In the case of ESO501, the simpl and cutoffpl models show good fits, with a χ² value of 422.19 for 416 degrees of freedom (reduced χ² ∼ 1.01). They show a local absorption column of N_H ∼ 0.07 × 10²² cm⁻², a cold disk temperature of ∼0.3 keV, a hot disk temperature of ∼1.1 keV, and a p index greater than 0.7, with an unabsorbed luminosity of 9 × 10³⁹ erg s⁻¹ for a distance of 7 Mpc in the 0.3–20 keV energy range.

For IC5052, both models also provide a good fit, with a χ² of 347.68 for 356 d.o.f. (reduced χ² ∼ 0.97), and yield with very similar spectral parameters (within uncertainties). Both models show a local absorption column of N_H ∼ 0.7 × 10²² cm⁻², a cold disk temperature close to 0.3 keV and a hot disk temperature around 1 keV with a p index of about 0.7. For both models, the unabsorbed luminosity, at a distance of 5.5 Mpc, is about ∼5 × 10³⁹ erg s⁻¹ in the energy range 0.5–20 keV.

In the spectrum of ESO501, both spectral models (simpl and cutoffpl) show an abrupt drop in emission for energies above 5 keV. In the simpl model, this drop is due to less efficient Compton scattering, indicated by a high spectral index (Γ = 4), resulting in a rapid reduction of high-energy photons. In the cutoffpl model, the decrease is due to the presence of a power-law tail with an exponential cutoff at 7.1 keV, which limits emission at higher energies.

For IC5052, both models (simpl and cutoffpl) provide a consistent interpretation of the spectrum, showing stable emission for energies above 3 keV. The simpl model, which incorporates Compton scattering of disk photons, achieves this with an index of Γ = 2.1 and a significant scattering fraction (c_f = 0.8), effectively extending the emission to higher energies. Similarly, the cutoffpl model, which includes an exponential cutoff at 7.1 keV, maintains a stable high-energy tail, consistent with the expected modulation of the emission at these energies. Despite their different physical foundations – Compton scattering in simpl and emission modulation in cutoffpl – both models fit the observed data well, reflecting the robustness of the thermal and high-energy components in IC5052, without indicating any significant differences in the underlying emission mechanisms.

3.3. Constrains on the nature of the compact object

In this section we perform a comparative analysis of the ULX sources ESO501 and IC5052 using several diagnostic tools. By examining the spectral characteristics of these sources in the context of the broader ULX population, we aim to constrain the nature of the compact objects driving the observed emission. Specifically, we analyze and discuss the position of these sources in a hardness-luminosity diagram (HLD), investigate their spectral properties using a color-color diagram (CCD), and determine the dominance of high-energy emission through Flux Ratio Analysis.

The results of our spectral and timing analysis indicate that both sources display characteristics consistent with super-Eddington accretion. Although no significant pulsations were detected in either source, with an upper limit for the PF below 11%, this absence of pulsations does not rule out the possibility that the systems contain NSs. For instance, both M82 X-2 and NGC 1313 X-2 have shown pulsations with a PF often below 11% (Bachetti et al. 2014; Sathyaprakash et al. 2019), suggesting that a low PF is not a necessary condition to exclude the presence of an NS in such systems.

3.3.1. Hardness-luminosity diagram analysis

In order to further explore the possible pulsating nature of these sources, we first constructed a HLD, similar to those used in studies of X-ray binaries and ULXs (see, e.g., Done & Gierliński 2003; Sutton et al. 2013; Gúrpide et al. 2021b). The HLD involves the calculation of fluxes and unabsorbed luminosities within a given energy range. We used the best fit parameters associated with the simpl model as described in Table 1 for ESO501 and IC5052.

We determined the total unabsorbed luminosity in the 0.3–10 keV energy range. We then computed the hardness ratio between the unabsorbed flux in a soft band (0.3–1.5 keV) with that in a hard band (1.5–10 keV). This distinction between soft and hard bands arises because the pulsating component in ULXs (PULXs) typically dominates at higher energies, while the soft component in ULXs tends to diminish around ∼1 keV (see Israel et al. 2017; Walton et al. 2018a).

According to the above definitions, ESO501 has a luminosity of $8.7_{-0.3}^{+0.4}\times {10}^{39}$ erg s⁻¹ and a hardness ratio of 1.5 ± 0.1. Although these values are similar to those observed in certain states of NGC 1313 X-1 (Sutton et al. 2013; Gúrpide et al. 2021a), as shown in Fig. 2, a direct comparison must be made with caution. While the two sources have comparable hardness ratios and luminosities in certain observations, their spectral properties are different. In particular, NGC 1313 X-1 exhibits a broader spectral variability, transitioning between hard and soft states (Sutton et al. 2013), while ESO501 has a more consistent spectral profile. In light of these considerations, although the position of ESO501 in the HLD suggests a similarity to the HUL state observed in NGC 1313 X-1, the spectral differences suggest that ESO501 does not fully correspond to a typical HUL state. The data instead indicate a possible BD state with contributions from both soft and hard spectral components. This interpretation is supported by the relatively high hardness ratio, which could result from a mixed contribution of thermal and nonthermal emission.

Fig. 2. HLD and CCD for ULXs taken from Gúrpide et al. (2021a) and Pintore et al. (2017), respectively. Left panel: Several ULX sources with different time periods, including PULXs and non-PULXs. All luminosities and fluxes are shown without absorption effects. The lightly shaded red and purple data points correspond to the sources studied in this work. The dashed black lines represent 10 and 100 times the Eddington limit for an NS with a canonical mass of 1.4 M⊙ (∼2 × 1038 erg s−1). Right panel: CCD generated from flux ratio calculations in the 2–4 keV, 4–6 keV, and 6–30 keV energy ranges using the optimal fits described in Table 1. The lightly shaded silver dots represent the sources analyzed in this study.

In the case of IC5052, the source has a luminosity of $4.2_{-0.4}^{+1.7}\times {10}^{39}$ erg s⁻¹ and a hardness ratio of $1.4_{-0.3}^{+1.2}$. While the nature of the compact object remains uncertain, it is located close to well-known objects such as Circinus ULX5, which could host either a BH or an NS (Mondal et al. 2021), and NGC 1313 X-2, a confirmed PULX. This leaves open the possibility that it could be a BH, a non-pulsating NS, or a PULX where the pulsations may intermittently appear and disappear. This transient nature means that pulsations are not consistently detectable throughout observations.

The results from the HLD analysis provide a preliminary classification of the accretion states of ESO501 and IC5052. To further explore the spectral properties of these sources, we next construct a CCD to examine how their high-energy components compare with those of other ULXs.

3.3.2. Color-color diagram analysis

To better understand the spectral differences in the high-energy components of these sources, we constructed a CCD. We derived the colors using fluxes in the 2–4 keV, 4–6 keV and 6–30 keV energy bands, obtained from the best fits shown in Table 1 for ESO501 and IC5052. This CCD is presented in Fig. 2, where hardness is defined as the ratio (6 − 30 keV)/(4 − 6 keV), and softness as (2 − 4 keV)/(4 − 6 keV).

In addition, we have included other sources from the work of Pintore et al. (2017) for comparison with the sources in our study, shaded in light gray. According to Pintore et al. (2017), sources classified as PULXs tend to occupy regions with harder spectra in the CCD, particularly when pulsations are detected, as seen in objects such as NGC 779 P13 and NGC 5907 X-1 (see Fig. 2). In contrast, ULXs likely to host BHs generally have softer spectra, although there is some overlap in the hardness ratios between different ULX types.

The position of ESO501 in the CCD (Fig. 2) places it in a region where both PULXs and ULXs without detectable pulsations are located. Although its position suggests that it is unlikely to be a PULX, the possibility of transient pulsations cannot be entirely ruled out, especially given that some sources, like NGC 1313 X-2, have shown pulsations only in certain observations. Furthermore, other sources with similar spectral properties, such as NGC 1313 X-1 and NGC 5408 X-1, have also not shown detectable pulsations, although the possibility that they host an NS cannot be ruled out. Therefore, while the spectral characteristics of ESO501 suggest it may host a non-pulsating NS or a BH, we cannot entirely exclude the possibility that it could exhibit transient pulsations.

In contrast, IC 5052, as shown in Fig. 2, is situated in a region of the CCD that is consistent both with known PULXs, such as NGC 7793 P13 and NGC 5907 X-1, and with other ULXs (e.g., Ho IX X-1, IC 342 X-1, NGC 1313 X-1) that have not shown detectable pulsations but are still considered potential candidates to host an NS with transient pulsations. This positioning supports the hypothesis that IC 5052 could be pulsating transiently.

Following the insights gained from the CCD analysis, we proceed to evaluate the dominance of high-energy emission through a flux ratio analysis to further constrain the nature of the compact objects in ESO501 and IC5052.

3.3.3. Flux ratio analysis

Building on the previous analyses, we calculate the ratio between the total flux in the 0.3–30 keV energy range and the flux modeled by cutoffpl in the same energy range. This ratio allows us to assess the dominance of high-energy emission–a characteristic feature often linked to PULXs and potentially indicative of a significant accretion column (Walton et al. 2018c).

For ESO501, the flux ratio is F_cutoffpl/F_tot ≲ 9%, indicating that the high-energy emission does not significantly dominate the total flux. This result suggests that the source may not be a PULX. Alternatively, it could be an NS in which pulsations are undetectable due to variations in the mass accretion rate, changes in system geometry, shifts in magnetic field configuration, among other possibilities (see, e.g., Bachetti et al. 2020; Imbrogno et al. 2024). In contrast, for IC5052, the flux ratio is F_cutoffpl/$F_{\mathtt{tot}}\sim {43}_{-9}^{+8}$%, suggesting a greater contribution from the high-energy emission. Although this value is below the typical PULX threshold (∼50%), could imply that IC5052 is also an NS with undetectable pulsations.

The comprehensive analysis in this section – including spectral characteristics, hardness ratios, and flux dominance – indicates distinct underlying mechanisms for ESO501 and IC5052. Assuming that these ULXs host BHs, we proceed to model their X-ray emission using a physical model that interprets thermal components within the framework of super-Eddington accretion modulated by system geometry and radiation-driven winds.

4. Physical interpretation

We now present a more physical model to account for the X-ray observations of both ULXs, based on that presented by Abaroa et al. (2023) and recently extended by Combi et al. (2024). Although the nature of the compact objects in the binary systems is unclear, in this model we assume that they are stellar-mass BHs. Under this assumption, we aim to constrain parameters such as the BH masses, accretion rates, and the photosphere of the super-Eddington disk-driven wind, and to characterize the electron population in the ULXs funnels.

4.1. Disk and wind models

We assumed for both ULXs that the BHs accrete matter at super-Eddington rates (M_input > Ṁ_Edd) from an evolved massive star. The Eddington rate depends only on the mass of the BH,

$$\begin{array}{c}\hfill M_{\mathrm{Edd}}=\frac{L_{\mathrm{Edd}}}{\eta c^2}\approx 2.2\times {10}^{-8}{M}_{\mathrm{BH}}{\mathrm{yr}}^{-1}=1.4\times {10}^{18}\frac{M_{\mathrm{BH}}}{M_{\odot }}\mathrm{g}{\mathrm{s}}^{-1},\end{array}$$(1)

with L_Edd the Eddington luminosity (defined as the luminosity required to balance the attractive gravitational pull of the BH by radiation pressure), η ≈ 0.1 the accretion efficiency, and c the speed of light.

The critical radius of the accretion disk is given by r_crit ∼ 400 ṁ_g, where ṁ = Ṁ_input/Ṁ_Edd and r_g = GM_BH/c² is the gravitational radius of the BH, with G the gravitational constant. This radius corresponds to the distance from the BH at which the outer standard disk (Shakura & Sunyaev 1973) changes to the inner disk dominated by radiation (Fukue 2004). The disk becomes geometrically thick in the inner region, where the ejection of winds by the radiation force helps regulate the mass accretion rate onto the BH (Ṁ_acc) at the Eddington rate.³ Such a regulation results in a total mass loss rate in winds Ṁ_w that is approximately equal to the accretion input, Ṁ_w ≈ Ṁ_input.

A self-similar solution for optically thick supercritical disks with mass loss can be found for the effective disk temperature, as in the appendix of Fukue (2004). The solution reads

$$\begin{array}{c}\hfill \sigma T_{\mathrm{eff}}^4=\{\begin{array}{c}\frac{3G{M}_{\mathrm{BH}}{\dot{M}}_{\mathrm{input}}}{8\pi r_{\mathrm{d}}^3}{f}_{\mathrm{in}},r_{\mathrm{d}}>r_{\mathrm{crit}}\hfill \\ \hfill \\ \frac{3}{4}\sqrt{c_3}\frac{L_{\mathrm{Edd}}}{4\pi r_{\mathrm{d}}^2},r_{\mathrm{d}}\le r_{\mathrm{crit}}\hfill \end{array},\end{array}$$(2)

with $\sqrt{c_3}=H/r_{\mathrm{d}}=tan\delta$, where H is the scale height of the disk, δ is the disk opening angle, and f_in = 1 − r_in/r_d ≈ 1 (since r_d > r_crit, then r_d ≫ r_in). Here, r_in is the location of the inner edge of the disk. The coefficient c₃ depends on the advection parameter, the adiabatic index of the gas, and the viscosity (see the appendix in Fukue 2004, for details). We note that in the critical region of the disk T_eff ∝ r^−1/2, contrary to what is expected for a standard disk (Watarai et al. 2000).

The effects of the radiation on the matter located at the escape surface of the disk in the region r_d ≤ r_crit have been calculated in detail by Abaroa et al. (2023). The radiation field transfers both energy and angular momentum to the wind, which escapes the system with a strong equatorial component, forming a funnel above the BH (Abaroa et al. 2023). This scenario is typical of ULXs under the assumption of super accretion (King 2009; Fabrika et al. 2021; King et al. 2023).

The wind itself is opaque until it reaches the photospheric radius, where it becomes transparent to its own radiation. This photosphere is defined by the condition that the optical depth, τ_w, is unity for an observer at infinity. If the velocity of the wind is relativistic, the optical depth in the observer’s frame depends in general on the magnitude of the velocity and the viewing angle. The position of the apparent photosphere from the equatorial plane, z_photo, is given by (Abaroa & Romero 2024)

$$\begin{array}{c}\hfill {\tau }_{\mathrm{w}}={\displaystyle {\int }_{z_{\mathrm{photo}}}^{\infty }}{\gamma }_{\mathrm{w}}(1-\beta cos\theta ){\kappa }_{\mathrm{co}}{\rho }_{\mathrm{co}}\mathrm{d}z=1,\end{array}$$(3)

where γ_w is the wind Lorentz factor, κ_co is the opacity in the comoving frame, and ρ_co is the comoving wind density. Since we assumed a fully ionized wind, the opacity is dominated by free electron scattering (κ_co = σ_T/m_p).

4.2. Disk emission

When the system is viewed face-on (i.e., with an inclination angle of i ≈ 0°), the only X-ray emission visible from the disk is that which escapes through the funnel (i.e., the emission produced near the innermost part). The radiation from the rest of the disk is absorbed and reprocessed in the wind and shifted to lower energies (see Abaroa et al. 2023, 2024). The observed energy-dependent luminosity of the disk is then L_observed(E_γ) = L_emitted(E_γ)⋅e^{−τ_w(E_γ)} (where e^{−τ_w(E_γ)} is the attenuation factor), so the disk spectrum is attenuated for energies lower than those corresponding to the radius at which the wind column above the disk has an optical depth τ_w ≳ 1.

The dense wind surrounds a low-density duct producing a geometrical beaming effect (Lasota & King 2023), which causes an observer to infer an isotropic bolometric disk luminosity given by

$$\begin{array}{c}\hfill L_{\mathrm{beamed}}\approx L_{\mathrm{Edd}}[1+\frac{ln\dot{m}}{b}],\end{array}$$(4)

where b is the beaming factor (King 2008). This factor varies with the normalized accretion rate as b ∝ ṁ⁻² for ṁ ≳ 8 (King et al. 2023). We note that all the emission from the disk reaching the observer is produced in its innermost region, as the dense wind absorbs the remaining radiation from the outer part. The beaming factor will then uniformly affect the disk spectrum on all energies.

4.3. Results

4.3.1. Thermal spectra

Figure 3 (left) shows the results for ESO501, where we plot the disk emission for a stellar-mass BH (solid blue line). Most of the data from this source are consistent with the radiation from the inner disk of a super-Eddington BH of 10 M_⊙ with an accretion 10Ṁ_Edd, following the model above, which is also consistent with the observational hypothesis of a ULX in a broadened state. This accretion rate leads to a beaming factor of b = 0.73. The opaque wind suppresses the disk emission for energies ≲300 eV. We also show for comparison the application of our model assuming alternative values for the accretion rates: 5Ṁ_Edd (dash-dotted lines) and 15Ṁ_Edd (dashed lines).

Fig. 3. Spectral energy distributions (SEDs) of the X-ray emission of ESO501 (left) and IC5052 (right) on a logarithmic scale. In both cases, red data correspond to the PN camera of XMM-Newton and purple data to the FPMA camera of NuSTAR. The solid blue lines correspond to the attenuated thermal emission from the inner accretion disk, and the orange line is the nonthermal Comptonization of the thermal photons by the relativistic electrons in the ULX funnels. The solid black line corresponds to the total emission (sum of thermal and nonthermal contributions) assuming the parameter values listed in Table 2. The disk emission of an IMBH is shown with a green line for IC5052. We also show, for comparison, the application of our model assuming alternative values for the accretion rates: 5ṀEdd (dash-dotted lines) and 15ṀEdd (dashed lines) for ESO501, and 5ṀEdd (dash-dotted lines) and 10ṀEdd (dashed lines) for IC5052.

The results for IC5052 are shown in Fig. 3 (right). The soft X-ray data are in agreement with a BH of 8 M_⊙ and an accretion of 7Ṁ_Edd (solid blue line), so beaming is not relevant in this case. The opaque wind suppresses the disk emission for energies ≲500 eV. As can be seen from the figure, the model for the sub-Eddington IMBH of 133 M_⊙ considered here is not in agreement with the data (solid green line). The more massive the BH, the cooler the innermost stable orbit of the accretion disk. On the other hand, a super-Eddington IMBH would surpass the X-ray luminosity observed. We also show for comparison the application of our model assuming alternative values for the accretion rates: 5Ṁ_Edd (dash-dotted lines) and 10Ṁ_Edd-pagination (dashed lines).

4.3.2. Nonthermal spectra

In both sources the emission extends to energies higher than ∼10 keV that cannot be of thermal origin from the disk; the hard X-rays found in our observations require an additional component. Such a component could be a hot plasma in the funnel above the BH (Combi et al. 2024). Since the radiation pressure in the funnel above the BH exceeds the gravitational pull, the plasma cannot be at rest and must be constantly removed and replenished. Since heating the gas to Comptonization temperatures requires time for thermalization, we prefer to assume a nonthermal component flowing through the funnel (Combi et al. 2024). This plasma can be heated by magnetic reconnection, and a nonthermal particle population could be created by turbulent acceleration (del Valle et al. 2016). The mildly relativistic electrons in the plasma then Comptonize the softer photons from the disk, producing the nonthermal tail. We modeled the Comptonization assuming in both cases an spectral index p = 2, a hadron-to-lepton ratio a = 0.1, and a fraction q_rel = 0.15 of the initial energy budget that goes to accelerate the relativistic particles. The solid orange line in Fig. 3 represents the Comptonization for these parameter values. The application of our model assuming alternative values for the accretion rates is shown in the figure in the same way as for the thermal components.

4.4. Final remarks

As can be seen from the spectral energy distributions in Fig. 3, the model is in good agreement for both ULXs with the available data for the set of parameters presented in Table 2. The main difference between the two sources from this work is that Comptonization plays a leading role in the case of IC5052, whereas in ESO501 the emission is strongly dominated by the thermal component. This supports the classification of the ULXs states as HUL and BD, respectively. Alternative values for the accretion rate considered in this paper underestimate (lower rates) and overestimate (higher rates) the data for both the thermal and nonthermal components.

The hypothesis of a sub-Eddington IMBH in IC5052 is hard to sustain, as can be seen in our results. The disk of a 133 M_⊙ BH has its emission peak at ∼1  keV and a luminosity of 10³⁷ erg s⁻¹ at 10 keV, while the observed power in X-rays at that energy is about two orders of magnitude higher. Furthermore, since there are no radiation-driven winds in a sub-Eddington system, the emission from the accretion disk is not absorbed, so in our model the IMBH luminosity overestimates the data for energies ≲3  keV.

Concerning the Comptonization, the conditions in the transparent funnel are not expected to be very stable, due to various instabilities, clumping of the wind, and magnetic reconnection (Combi et al. 2024). Hence, some degree of rapid variability could be present. Longer observations will test this in real time.

5. Conclusions

In this study we have thoroughly analyzed the ULX sources ESO501 and IC5052 using combined observations from XMM-Newton and NuSTAR. The results of our spectral and timing analysis suggest that both sources exhibit properties consistent with super-Eddington accretion. Although no significant pulsations were detected in either source, with a pulsed fraction below 11%, we cannot rule out the possibility that they host NSs with transient pulsations.

The spectral analysis indicates that ESO501 could be in a BD state, while IC5052 appears to be in an HUL state. The HLD and CCD support this interpretation, although the higher-energy emission of IC5052 could indicate a transition state or a complex accretion geometry. The flux ratio analysis shows that high-energy emission is not dominant in ESO501, which suggests the absence of an accretion column. In IC5052, however, the higher flux ratio implies a contribution from the accretion column, possibly due to a magnetic field associated with an NS. Nonetheless, the absence of pulsations supports a scenario in which the host is either a BH or an NS with transient pulsations. To interpret the thermal and nonthermal emission components observed in these sources, under the assumption that the compact object is a BH, we also applied a physical model developed for super-Eddington accretion.

Confirming the true nature of these compact objects will require additional observations, particularly those sensitive to pulsations, as pulsation detection remains the only conclusive method for identifying a PULX. Next-generation telescopes with enhanced observational capabilities, such as Athena, will be crucial for advancing our understanding of these ULX sources.

Acknowledgments

Funded by the European Union (Project 101183150 – OCEANS). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Executive Agency (REA). Neither the European Union nor REA can be held responsible for them. FAF, JAC, and FG acknowledge support by PIP 0113 (CONICET). EAS acknowledges support by the Spanish Agencia estatal de investigación via PID2021-124879NB-I00. FAF is a postdoc fellow of CONICET. JAC, GER, and FG are CONICET researchers. FG acknowledges support from PIBAA 1275 (CONICET). JAC was supported by Consejería de Economía, Innovación, Ciencia y Empleo of Junta de Andalucía as research group FQM-322. FG, JAC and GER were also supported by grant PID2022-136828NB-C42 funded by the Spanish MCIN/AEI/ 10.13039/501100011033 and “ERDF A way of making Europe”. GER acknowledges financial support from the State Agency for Research of the Spanish Ministry of Science and Innovation under grants PID2019-105510GB-C31AEI/10.13039/501100011033/, and through the “Unit of Excellence María de Maeztu 2020-2023” award to the Institute of Cosmos Sciences (CEX2019-000918-M). Additional support came from PIP 0554 (CONICET).

References

All Tables

All Figures

	Fig. 1. XMM-Newton+NuSTAR unfolded spectra and residuals of ESO501 (left column) and IC5052 (right column). Each model and best-fit statistic are specified.
In the text

Fig. 2. HLD and CCD for ULXs taken from Gúrpide et al. (2021a) and Pintore et al. (2017), respectively. Left panel: Several ULX sources with different time periods, including PULXs and non-PULXs. All luminosities and fluxes are shown without absorption effects. The lightly shaded red and purple data points correspond to the sources studied in this work. The dashed black lines represent 10 and 100 times the Eddington limit for an NS with a canonical mass of 1.4 M⊙ (∼2 × 1038 erg s−1). Right panel: CCD generated from flux ratio calculations in the 2–4 keV, 4–6 keV, and 6–30 keV energy ranges using the optimal fits described in Table 1. The lightly shaded silver dots represent the sources analyzed in this study.

In the text

Fig. 3. Spectral energy distributions (SEDs) of the X-ray emission of ESO501 (left) and IC5052 (right) on a logarithmic scale. In both cases, red data correspond to the PN camera of XMM-Newton and purple data to the FPMA camera of NuSTAR. The solid blue lines correspond to the attenuated thermal emission from the inner accretion disk, and the orange line is the nonthermal Comptonization of the thermal photons by the relativistic electrons in the ULX funnels. The solid black line corresponds to the total emission (sum of thermal and nonthermal contributions) assuming the parameter values listed in Table 2. The disk emission of an IMBH is shown with a green line for IC5052. We also show, for comparison, the application of our model assuming alternative values for the accretion rates: 5ṀEdd (dash-dotted lines) and 15ṀEdd (dashed lines) for ESO501, and 5ṀEdd (dash-dotted lines) and 10ṀEdd (dashed lines) for IC5052.

In the text