The high energy X-ray probe (HEX-P): studying extreme accretion with ultraluminous X-ray sources

1 Introduction

Ultraluminous X-ray sources (ULXs, see Kaaret et al., 2017; Fabrika et al., 2021; King et al., 2023; Pinto and Walton, 2023 for recent reviews) are off-nuclear X-ray sources whose apparent luminosities exceed the Eddington limit for a stellar-mass black hole (e.g., $\sim 10^{39}$ erg s⁻¹ for a $\sim 10M_{\odot }$ black hole).

ULXs commonly have X-ray spectra consisting of two thermal components; the lower-energy component has a characteristic temperature of $<1$ keV, whereas the higher-energy component shows a cut off at $>3$ keV. The low-energy component was initially interpreted as evidence for standard, sub-Eddington accretion onto intermediate-mass black holes (Miller et al., 2003). However, the cutoff, hinted at in high-quality observations with XMM-Newton since the mid-2010s (Roberts et al., 2005; Stobbart et al., 2006; Middleton et al., 2015a), was interpreted as Comptonization of disk photons by an optically thick corona and super Eddington accretion (e.g., Gladstone et al., 2009). Hard X-ray coverage provided by NuSTAR (Harrison et al., 2013) gave the highest-significance detections of this higher-energy cutoff, confirming the likely super-Eddington nature of most ULXs (e.g., Bachetti et al., 2013; Walton et al., 2013; Rana et al., 2015). Following the first detection of pulsations in M82 X-2 with NuSTAR (Bachetti et al., 2014), a growing number of pulsating ULXs (PULXs, also referred to as ultraluminous pulsars, ULPs; e.g., Fürst et al., 2016; Israel et al., 2017; Carpano et al., 2018; Rodríguez Castillo et al., 2020) have been found, showing that neutron stars (NSs) are able to radiate at hundreds of times their Eddington limit. The pulsations tend to be associated with an additional hard spectral component above 10 keV, which is reminiscent of the hard, curved spectra found in accreting pulsars (e.g., Pintore et al., 2017; Walton et al., 2018b).

The detection of some Galactic super-Eddington NSs provide a possible link between high-luminosity X-ray binaries (XRBs) and ULXs (Wilson-Hodge et al., 2018, e.g., Swift J0243+6124); in general, studying the highest-luminosity end of Galactic XRBs in the Milky Way Galaxy and external galaxies might provide clues on the onset of the ULX regime [see Fabbiano (2006) for a review, and Connors et al. (2023), in this volume, for HEX-P’s contribution to studies of this Galactic population].

Except for the few cases where pulsations have been observed, implying a neutron star (NS) accretor, determining the nature of the compact object in ULXs remains a major challenge for current observatories and a key open question for the vast majority of sources. One way to address this issue is through population synthesis studies, which suggest that NS-ULXs may dominate the ULX population (King et al., 2017; Middleton and King, 2017; Wiktorowicz et al., 2019; Khan et al., 2021). PULXs can reach luminosities above 10⁴¹ erg s⁻¹ (e.g., NGC 5907 X-1, Israel et al., 2017), and given the similarity of PULX spectra to much of the ULX population, NSs might indeed power the majority of ULXs (Pintore et al., 2017; Walton et al., 2018b). In addition, a very small number of ULXs have shown evidence for cyclotron resonance scattering features (CRSFs) in their spectra¹ (e.g., Walton et al., 2018a; Brightman et al., 2018), providing another route to identifying the presence of a NS, as well as a direct estimate of the magnetic field strength (likely indicating a weaker dipole but stronger multipolar field) close to the NS surface (Middleton et al., 2019a; Kong et al., 2022).

The mechanism for the very high luminosities of ULXs is still debated. It is likely that the emission is partially collimated by an optically thick, radiatively driven outflow (King et al., 2001) launched from the large scale-height disk expected at very high mass-transfer rates (Shakura and Sunyaev, 1973). Indeed, many ULXs show evidence of high-velocity winds and outflows likely inflating the $\sim 100$ pc interstellar bubbles found around many ULXs at various wavelengths (Stobbart et al., 2006; Middleton et al., 2014; 2015b; Pinto et al., 2016; Pakull and Mirioni, 2002; Gúrpide et al., 2022; Belfiore et al., 2020, see Section 2.1). However, these mass-loaded winds likely affect our ability to locate pulsing ULXs (e.g., Mushtukov et al., 2017), and collimation alone is difficult to reconcile with the high pulsed fraction concurrent with the high luminosities of PULXs (e.g., Israel et al., 2017). Another possibility is that a locally strong, possibly non-dipolar, magnetic field is capable of altering the local Thomson cross section (Basko and Sunyaev, 1976; Dall’Osso et al., 2015; Ekşi et al., 2015; Mushtukov et al., 2015; Israel et al., 2017; Brice et al., 2021). These processes might also operate together to a greater or lesser extent.

Most ULXs are known to have other bright X-ray sources within 1′, in particular those located outside the Local Group (see Section 3.1.2). A clear detection and isolation of sources is therefore not possible given the low resolution power of the high energy detectors in use, impeding a profound and detailed study of the whole population of ULXs. The understanding of ULXs could reach considerable advancement by the combination, in one instrument, of high angular resolution, high time resolution and sensitivity in the hard X-rays.

The High-Energy X-ray Probe (HEX-P; Madsen et al., 2023) is a probe-class mission concept that offers sensitive broad-band X-ray coverage (0.2–80 keV) with an exceptional combination of spectral, timing and angular resolution capabilities. It features two high-energy telescopes (HETs) that focus hard X-rays, and one low-energy telescope (LET) that focuses lower energy X-rays.

The LET consists of a segmented mirror assembly, coated with iridium on monocrystalline silicon that achieves a half power diameter (HPD) of 3.5″, and a low-energy DEPFET detector, of the same type as the Wide Field Imager (WFI; Meidinger et al., 2020) onboard Athena (Nandra et al., 2013). It has 512 × 512 pixels that cover a field of view of 11.3′ × 11.3′, an effective passband of 0.2–25 keV, and a full frame readout time of 2 ms, which can be operated in a 128 and 64 channel window mode for higher count-rates, to mitigate pile-up and allow for faster readout. Pile-up effects remain below an acceptable limit of $\sim 1\%$ for fluxes up to $\sim 100$ mCrab–a lot higher than typical extragalactic ULX fluxes–in the smallest window configuration. The high angular resolution and low background result in a factor ∼2 sensitivity improvement with respect to XMM-Newton.

The HET consists of two co-aligned telescopes and detector modules. The optics are made of Ni-electroformed full shell, mirror substrates, leveraging the heritage of XMM-Newton (Jansen et al., 2001), and coated with Pt/C and W/Si multilayers for an effective passband of 2–80 keV. The high-energy detectors are of the same type flown on NuSTAR (Harrison et al., 2013), consisting of 16 CZT sensors per focal plane, tiled 4 × 4, for a total of 128 × 128 pixels spanning a field of view of 13.4′ × 13.4′. This improvement in angular resolution allows for a much more sensitive instrument configuration compared to NuSTAR.

The broad HEX-P X-ray passband and superior sensitivity compared with existing facilities provides a powerful and important opportunity to study ULXs across a wide range of energies, luminosities, and timescales.

2 Challenges and open questions

A number of major questions in the field of ULXs and super-Eddington accretion remain unanswered, preventing us from achieving a comprehensive understanding of their phenomenology and detailed accretion physics. We briefly describe them here.

2.1 Where does the transferred mass go?

Understanding the fraction of matter being expelled in winds versus the fraction accreted is key to estimate the possible evolutionary paths of massive binary systems–for instance, the formation of double compact object binaries and systems that create large ionising flux at high redshifts (Fragos et al., 2013)—and for understanding potential supermassive black hole (SMBH) growth rates in the early universe. A possible means for measuring the current mass loss rate is via the absorption lines, which are indicators of quasi-relativistic winds [see Pinto et al. (2016)], requiring observations at $<2$ keV. Many ULXs are associated with optical or radio bubbles (Pakull and Mirioni, 2002; Soria et al., 2021), which are sometimes as large as 400 pc in diameter and mechanically inflated (Gúrpide et al., 2022), allowing an estimate for the integrated kinetic energy released by the winds. The mass inflow rate can be crudely estimated from soft X-ray spectra [see Poutanen et al. (2007)] or, should the mass transfer dominate the angular momentum distribution within the system, from the change in the orbital period [e.g., see Bachetti et al. (2022), who estimate the mass transfer rate from a decrease in the orbital period of M82 X-2 through pulsar timing].

Winds in XRBs are commonly identified via absorption lines in the Fe K band from hot plasmas (typically Fe XXV-XXVI), and in the soft band from a forest of H- and He-like ions of several elements from C to Fe [for a review see Neilsen and Degenaar (2023)]. Initial searches for Fe XXV-XXVI in ULXs did not result in detections [see, e.g., Walton et al. (2013)], most likely due to their soft spectra or high ionisation state of the gas [for a review on ULX winds see Pinto and Walton (2023)]. Evidence of atomic lines in low-resolution ULX X-ray spectra were found in the form of residuals around 1 keV with respect to the continuum model through high-count-rate XMM-Newton EPIC spectra (Stobbart et al., 2006; Middleton et al., 2014; 2015b). Later on, high-spatial-resolution instruments onboard Chandra showed that the residuals were associated with the ULX itself rather than the interstellar medium in the host galaxy (Sutton et al., 2015), which was supported by the anticorrelation between their strength and the ULX X-ray spectral hardness (Middleton et al., 2015b). The first unambiguous proof of winds was obtained through long observations (300–500 ks) with the high-spectral-resolution reflection grating spectrometer (RGS) onboard XMM-Newton. In particular, Pinto et al. (2016) discovered a forest of emission lines (at rest-frame laboratory wavelengths) and absorption lines (mainly O VII-VIII and Ne IX-X), blueshifted by $\sim 0.2c$ in the two well-known ULXs NGC 1313 X-1 and NGC 5408 X-1, in excellent agreement with the predictions of winds driven by strong radiation pressure in super-Eddington accretion disks (Ohsuga and Mineshige, 2011; Jiang et al., 2014). Modelling of these lines permits an estimate of the kinetic power within the wind (affected by the uncertainties in the covering fraction), which matches well with the luminosities required to inflate the supersonically expanding ULX bubble nebulae [detected primarily in the optical and showing broad lines, see, e.g., recent work with the MUSE camera on the Very Large Telescope by Gúrpide et al. (2022)].

Features indicating a wind have been found in the majority of ULX RGS spectra with ≳ 1, 000 net counts (Pinto et al., 2017; Kosec et al., 2018; Pinto et al., 2021), as well as in the Chandra gratings observation of the transient Galactic super-Eddington NS, Swift J0243.6+6124 (van den Eijnden et al., 2019). The luminosity of the plasma producing the emission lines varies slightly over time and is very large (L_X ∼ 10³⁸ erg s⁻¹; Pinto et al., 2020a), orders of magnitude higher than in Eddington-limited Galactic XRBs (Psaradaki et al., 2018). Absorption lines at 8.6–8.8 keV from the hotter Fe K component of the wind were finally found (albeit at lower significance) in a few ULXs with hard spectra (Walton et al., 2016b; Brightman et al., 2022). The emission and absorption lines in the soft band (0.3–2 keV) have equivalent widths similar to those of Galactic XRBs [EW ∼ 5–10 eV, see Kosec et al. (2021)]. The Fe K absorption features, albeit rarer, appear deeper (EW ∼ 50–200 eV) and similar to the most extreme wind phases observed in Galactic XRBs (which have velocities of $\sim 0.05c$; e.g., King et al., 2012).

The very low detection rate of Fe K lines arises as a combination of low effective area in current instruments, high background, and low ionic fractions. The latter is a consequence of the soft radiation field of ULXs (e.g., Pinto and Walton, 2023). However, in ULXs likely seen at low inclinations [e.g., the hard ultraluminous, or HUL, systems of Sutton et al. (2013) with spectral slopes Γ < 2], the continuum is strong, likely overwhelming the lines in the soft band (e.g., Middleton et al., 2015b). The Fe K features might then be detected, as the irradiating continuum is hard, as long as the gas is not overionised in the wind cone and the instrumental background is sufficiently low at 7–9 keV.

Although absorption lines have proven elusive in HUL sources (such as Holmberg IX X-1 and most pulsating ULXs), should such features be detected, we could start to understand a) the 3D structure of the wind, and b) the actual mass outflow rate (and true kinetic power). The latter has important cosmological consequences for the net growth rate of SMBHs (e.g., Volonteri and Rees, 2005; Volonteri et al., 2015) and the overall role which ULXs play in their host galaxies, particularly at the peak of star formation (e.g., Prestwich et al., 2015). HEX-P will dramatically improve our ability to search for the Fe K wind component in ULXs thanks to its broader energy band (needed for modelling the continuum, which is crucial to search for shallow spectral features, see Figure 8), improved sensitivity compared to XMM-Newton and NuSTAR (for maximizing total source counts), and its narrower PSF (since ULXs tend to reside in crowded regions).

2.2 What is the content of the ULX population?

One of the open questions regarding the nature of ULXs is the prevalence of NS accretors and the strength of their magnetic fields. The detection of fast (≳ 1 Hz) pulsations unequivocally identifies a ULX as a NS, and pulsar timing can be used to indirectly estimate the magnetic field (e.g., Ghosh and Lamb, 1979; Wang, 1996). The detection of cyclotron lines in X-ray spectra allows us to tackle these two problems simultaneously, as it enables us to identify the presence of a NS where pulsations may be undetectable or absent (Brightman et al., 2018), while simultaneously offering an indication of the magnetic field strength (Walton et al., 2018a; Middleton et al., 2019a). To date, only two potential CRSFs have been reported in extragalactic ULXs (Walton et al., 2018a; Brightman et al., 2018), and one Galactic super-Eddington source, Swift J0243+6124 (Kong et al., 2022). The broadband coverage of HEX-P–and its higher sensitivity above 10 keV (where these lines are typically detected in NS XRBs; Staubert et al., 2019)—will allow us to increase the population of identified NS-ULXs by detecting CRSFs (see Section 3.3.3).

At the moment, due to the limited amount of high quality data, the sample of well-studied ULXs is very small (a few tens, compared to the thousands known; Walton et al., 2022; Tranin et al., 2023), restricting our ability to perform deep searches for pulsations and CRSFs, and hence preventing us from assessing the fraction of NSs in the ULX population. Given the energy dependence of the pulsing component (e.g., Brightman et al., 2016) and the potential for the CRSF to be located at a wide range of energies (depending on the magnetic field strength), broad-band observations, especially at high energies, are crucial.

A similar energy dependence also characterizes other forms of variability in addition to pulsations, such as quasi-periodic oscillations (QPOs), which can in principle be used to classify accreting objects (Lewin et al., 1988; Belloni and Hasinger, 1990; van der Klis, 2005). ULXs are known to show QPOs (Atapin et al., 2019), some of which match the phenomenology observed in other accreting systems, e.g., QPOs that have been detected over a broad-band [like M82 X-1’s 50–300 mHz QPO, see, e.g., El Byad (2021)] are more significant at higher energies (≳ 10 keV). The broadband spectro-temporal study of a rich sample of ULXs will help single out possible candidate black hole ULXs (e.g., Wolter et al., 2018). In particular, it will be crucial for studying the brightest ULXs, also referred to as hyperluminous X-ray sources (HLXs; Farrell et al., 2009; MacKenzie et al., 2023, still considered among the best candidates to host intermediate-mass black holes (IMBHs) although are rare and mostly found in very distant galaxies.

2.3 What drives the changes in ULX brightness/spectral shape?

ULXs have been observed to change dramatically in brightness on timescales of days, weeks, and months (e.g., Walton et al., 2016a; Gúrpide et al., 2021a; b). Possible causes of this variability include precession of the disk and wind (e.g., Pasham and Strohmayer, 2013; Middleton et al., 2015b; Luangtip et al., 2016; Amato et al., 2023) driven by a variety of processes [see the discussion in Middleton et al. (2018), (2019b)], changes in accretion rate [for instance by the launching of a thermal wind at large radius (Middleton et al., 2022), or by modulating the mass transfer rate in a radiatively driven stellar wind], or the onset of a propeller state associated with a NS close to spin equilibrium. The latter has been extensively searched for (Earnshaw et al., 2018; Song et al., 2020), with the strongest evidence seen in the changing spin period of NGC 5907 ULX-1 (Fürst et al., 2023). With a well considered monitoring strategy, it should also be possible to identify propeller states accompanying smaller changes in luminosity by isolating changes in the spectrum at high energies [see Middleton et al. (2023)], provided that a broad energy bandpass is available. Determining the mechanism behind the long timescale changes in ULXs then requires several key ingredients: a) high energy coverage and sensitivity to locate and track pulsations, b) low energy coverage and sensitivity to explore how the wind responds to changes in the broad band spectrum (e.g., Middleton et al., 2015b; Pinto et al., 2020b), and c) broad simultaneous coverage at high and low energies in order to locate the hard component switching off (whilst the soft emission remains stable), which is an indicator of a likely propeller transition, and more generally to detect spectral transitions (Amato et al., 2023).

3 Simulations

To demonstrate the improved capabilities of HEX-P compared to the current primary instruments for ULX studies, XMM-Newton and NuSTAR, we performed a detailed HEX-P simulation of NGC 253, a Galaxy with a large population of bright XRBs and ULXs, using the SImulation of X-ray TElescopes (SIXTE) software package (Dauser et al., 2019), with the procedure described in Section 3.1.1. This simulation provides precise estimates of the source and sky count rates, and the instrumental background, to set up reliable estimates for the detection and temporal and spectral characterization of sources in crowded fields. These estimates were refined through follow-up simulations, using the count rates and background from the SIXTE simulations as inputs.

All the simulations presented in this paper were produced with a set of response files that represent the observatory performance based on current best estimates as of Spring 2023 (see Madsen et al., 2023). The effective area is derived from a ray-trace of the mirror design including obscuration by all known structures. The detector responses are based on simulations performed by the respective hardware groups, with an optical blocking filter for the LET and a Be window and thermal insulation for the HET. The LET background was derived from a GEANT4 simulation (Eraerds et al., 2021) of the WFI instrument, and the HET background was derived from a GEANT4 simulation of the NuSTAR instrument; both simulations assume an L1 orbit for HEX-P. XMM-Newton and NuSTAR simulations were also performed for comparison, using the official set of response files distributed with SIXTE².

3.1 End-to-end simulations of an NGC 253-like galaxy

The Chandra Source Catalogue v.2.0 (Evans et al., 2010; 2020) contains 1583 ULXs [following the criteria from Walton et al. (2022)]. Among these ULXs, 14% have other X-ray sources within 5″, 45% within 20″, and 64% within 40′′. The vast majority of the ULXs appear close to the nuclear region of their host galaxy, particularly for the ULXs residing in more distant galaxies. In these conditions, improvements in the angular resolution, whilst not degrading the effective area, timing capabilities, and hard X-ray response, are key for studying ULXs.

ULXs are mostly extragalactic objects; some galaxies, such as M51 (7.2 Mpc; Terashima and Wilson, 2004), NGC 253 (3.5 Mpc; Wik et al., 2014), and the Cartwheel Galaxy (150 Mpc; Wolter et al., 2018), are known to contain a relatively large number of bright accreting sources and ULXs, and so far, only Chandra is able to fully resolve their positions. It is possible that a population of elusive, highly absorbed ULXs (e.g., Luangtip et al., 2015; West et al., 2023), has still escaped detection (or identification as ULXs), awaiting high-energy response (and an ability to spatially separate them).

With its high sensitivity and angular resolution, HEX-P will be able to detect and characterize the spectra of many new ULXs. Simulations indicate that the effective area and background of HEX-P will allow us to not only identify ULXs from their characteristic spectra, but also detect important spectral features such as emission and absorption lines associated with strong winds, or cyclotron resonance features that allow for an estimate of NS magnetic fields (see Figure 8). Moreover, these capabilities are key for detecting pulsations from pulsars in crowded environments. We demonstrate this in the Sections below.

3.1.1 SIXTE setup

For complete, end-to-end simulations, we used the SIXTE software package (Dauser et al., 2019), a Monte Carlo simulation toolkit for X-ray astronomical instrumentation. This software is able to take into account very accurately the source properties, including simulating complicated fields-of-view with diffuse emission, and characteristics of the optics and detection system of the mission. The SIXTE team worked with the leads of the HEX-P mission to provide up-to-date instrument specification files for HEX-P. The SIXTE website also contains hardware specifications for a number of existing and future missions that allow an easy comparison between them.

To run a simulation with SIXTE, we first generate a source input (SIMPUT-format) file containing: 1) an image of the diffuse emission; 2) its spectral shape defined as a Xspec-format xcm file; and 3) a list of point sources, having different spectral shapes and timing properties, including pulsations and/or aperiodic variability. Through the script run_sixt, we then generate realistic simulations for HEX-P, NuSTAR, and XMM-Newton.

3.1.2 Galaxy and source parameters

We set up a SIMPUT file using a map of the diffuse emission of NGC 253, and 70 point sources³.

The map of diffuse emission was obtained by removing all point sources from a deep Chandra map of the field. We assigned to it a thermal plasma spectrum [based on Wik et al. (2014); see the details in Lehmer et al. (2023), in this volume]. The total 0.5–7 keV flux of the diffuse emission was 9.6 ⋅ 10^–12 erg s⁻¹ cm⁻². The 70 point sources selected were the 70 brightest sources from the Chandra catalogue. We only conserved the position and flux between 0.5 and 8 keV, but otherwise we changed the spectral and temporal properties of each source, as follows:

• each source had a ULX-like cutoff-powerlaw spectrum, with random parameters uniformly distributed in intervals deemed plausible by comparison with the spectral characteristics of known ULXs (e.g., Pintore et al., 2017): cutoff E_cut between 1 and 10 keV, and Γ between −2 and 0;

• pulsations with a pulsed fraction between 20% and 100% [NGC 300 ULX has ∼80% pulsed fraction over the full band, see e.g., Carpano et al. (2018), Vasilopoulos et al. (2019)], a pulse period between 1 ms and 100 s, and a pulse profile described by a Von Mises distribution with κ between 1 and 5.

We then ran full SIXTE simulations for HEX-P, NuSTAR, and XMM-Newton to produce images such as those shown in Figure 1. To analyze the simulations and estimate count rates from the sources and various components of the background, we used circular source regions centered at the source position with a radius corresponding to half of the half-power diameter of the point spread function used by SIXTE for the simulation.

FIGURE 1. 1-Ms Simulation of a NGC 253-like starburst galaxy showing the various components including diffuse hot gas and XRBs. Here we compare the performance of XMM-Newton and NuSTAR [compare to Pietsch et al. (2001), Wik et al. (2014)] with the instruments onboard HEX-P, using the procedure described in Section 3.1.2. HEX-P is able to resolve most of the bright X-ray sources in nearby galaxies, enabling a detailed study of both the source and diffuse emission. In particular, HEX-P is able to resolve and study in detail all ULXs, including the ones with known notable features (pulsations, lines, eclipses; see text), with better sensitivity and equivalent-or-better energy and timing resolution compared to XMM-Newton and NuSTAR.

For each source region, we used the internal description of data in SIXTE to estimate how many photons came from the source n (SRC_ID = n), how many from the instrumental background (SRC_ID == −1) and how many from the sky background [(SRC_ID != −1)&(SRC_ID != n)].

3.2 Pulsar simulations

Assuming that one has averaged enough periodograms to yield normally-distributed powers, Lewin et al. (1988) estimates the significance of the detection of a broad feature against a Poisson background (e.g., a high frequency QPO) in a periodogram as

where I is the total count rate (sky, instrument, source), r the rms fractional variability of the QPO over the total flux, T_obs the observing time, and Δν the width of the feature. This means that the observing time necessary to detect a QPO scales with the inverse fourth power of the rms.

Pulsar searches are usually performed with single periodograms, or proxies of them such as the $Z_n^2$ statistic (Buccheri et al., 1983) or the Epoch Folding search (Leahy et al., 1983), for which one cannot assume normality of powers. However, it can be shown that a similar dependence on the rms is also valid in this case (Figure 2). Therefore, in order to perform a sensitive search for pulsations, the key ingredient is reducing all components of the background so that the rms amplitude is reduced to the intrinsic source rms. When pulsars are not dominating their fields-of-view, raw sensitivity of an instrument is less important than its ability to separate the target from the other sources. ULXs are borderline in this respect, as they are often relatively bright sources in their host galaxies, but being extragalactic, it is common to observe other ULXs or an active galactic nucleus (AGN) separated by less than a few arcminutes. The vast majority of known ULXs in nearby galaxies can be confidently separated by HEX-P (see Section 3.1). In addition, their hard X-ray tail, which is where the pulsed emission is strongest, is covered by the energy range where HEX-P is most sensitive. With the typical parameters of ULXs, HEX-P can detect pulsations in a fraction of the time it takes NuSTAR or XMM-Newton, and HEX-P can even outperform timing-dedicated, non-focusing instruments, which may have much larger effective areas but far poorer angular resolution.

FIGURE 2. Comparison of the single-trial detection significance (in number of sigma) of a quasi-sinusoidal pulsed profile with a $Z_n^2$ search, given a total number of counts and the pulsed flux fraction (expressed as the counts in the pulse profile divided by the total counts). The sensitivity degrades rapidly when the pulsar is not dominating the emission, and a factor of $\sim 5$ lower in pulsed fraction can make a pulsar undetectable. Dashed lines show the expected dependence from a law that scales like Eq. 1.

3.2.1 Simulation setup

To simulate time-variable fluxes, we use inverse transform sampling. By appropriate normalization, a (positive definite) time series can be considered a probability density function. We first calculate its cumulative distribution function (CDF) by simple integration. Then, we generate random values of the CDF by extracting a series of values uniformly distributed between 0 and 1. We then invert this CDF by calculating (through any interpolation algorithm) the time values corresponding to that value of the CDF. This can be applied to a variable aperiodic light curve, for example, generated through the Timmer and Koenig (1995) method, or to the phases of a pulsar given its pulse profile. If we start from pulsar phases, further corrections to these phases can be made to accommodate spin derivatives, with the usual Taylor series expansion $\mathrm{\Delta }\varphi =\dot{\nu }t+1/2\ddot{\nu }{t}^2+\dots$. We can further generate integer pulse number values, sum them to the pulse phases, and transform them into times by multiplying by the pulse period. These methods are already implemented in the stingray (Huppenkothen et al., 2019) and HENDRICS (Bachetti, 2018) packages that we use for standalone timing simulations, while the SIXTE package takes care of end-to-end simulations (Dauser et al., 2019, see below).

3.2.2 A note on frame time

The frame time of CCD-like detectors reduces the sensitivity to pulsations with periods below 10–50 frame times depending on the profile shape (because of the spread of each frame over multiple phases of the pulsations). Hence, we can expect that the 73 ms frame time in full window of the EPIC-pn camera onboard XMM-Newton will greatly reduce sensitivity to sinusoidal pulsations with periods ≲ 500 ms, and sharp pulsations (e.g., some rotation-powered pulsars) with even much longer periods. In Figure 3, we generated ∼100,000 photons with phases following a Von Mises distribution of concentration (κ = 10), transformed the phases into arrival times by adding an integer pulse number and multiplying by the period, and applied a fixed frame time to the event times. Then, we folded the profiles using the standard procedure of adding a random number uniformly distributed around ±Δt/2 to avoid artifacts⁴. The profiles are clearly distorted, losing sharpness and height as the pulse period approaches the sampling period, until the pulsation is not detectable anymore.

FIGURE 3. Distortion of a pulse profile when the pulse period p is comparable to the frame time Δt. This is relevant for CCD-based detectors like XMM-Newton’s cameras (5.7 ms frame time in small window mode) or HEX-P/LET (2 ms frame time).

XMM-Newton has special data acquisition modes that reduce the frame time while sacrificing field-of-view. For example, the “small window” mode reaches a $\sim 5.7$ ms frame time by using only $\sim 1/40$th of the pixels. Moreover, the “timing” (0.03 ms) and “burst” (7 μs) modes also sacrifice angular resolution in order to further increase timing precision (Strüder et al., 2001). In contrast, the LET onboard HEX-P has a 2 ms frame time in its standard full frame operation mode, which allows the detection of significantly shorter pulse periods than XMM-Newton, even when conducting surveys or in crowded fields-of-view, making HEX-P a game-changer for extragalactic pulsar searches. The CdTe/CZT-based high-energy instruments such as NuSTAR and the HET onboard HEX-P have a $\sim 10\mu$s timing resolution, which is more than enough to detect all known periodic and aperiodic phenomena expected in compact objects (≲2,000 Hz). In the following, we assume that all pulsars have periods detectable by all the instruments that are being compared (≳1 s).

3.2.3 Pulsar searches

We carried out two different procedures to search for (and characterize) pulsations in the simulated data described in Section 3.1, one blind search (to evaluate possible confusion of different unresolved candidates) and one directed at known pulsations (using the input periods, and evaluating the increase of detection significance with observing time).

For blind pulsation searches, we selected events from all the known source positions, using extraction regions equal to the HPD of each instrument. We then ran a standard Fourier search with a periodogram, to obtain a list of candidate frequencies (Lovelace and Sutton, 1969). Depending on the source position and the instrument used, there was significant contamination between the sources; this made frequency candidates from nearby sources appear in the search of others. This of course does not represent a problem. Around each pulsar candidate, we ran zoomed searches using the quasi-fast folding algorithm (Bachetti et al., 2020) as implemented in HENDRICS (Bachetti, 2018). This method allows a fast exploration of the frequency-frequency derivative plane around pulse candidates. As a metric for pulse detection we used the H test (de Jager et al., 1989) modified for binned profiles (Bachetti et al., 2021). We are mostly interested in the comparison between the instruments, and all parameters affecting the statistics are assumed the same for each⁵ (principally the observing time, which increases the frequency resolution of the FFT and hence the number of trials). We selected a H test statistic of 50 in any instrument as criterion to declare a rough detection. This resulted in $\sim 16$ out of 70 sources being detected by one instrument or the other, and in particular some of the ones in the very center of the galaxy. For all these detections, we plot the value of the H test calculated from all instruments at the best frequency from the search in Figure 4. The results show that HEX-P’s instruments provide a comparable or higher sensitivity to pulsations compared to XMM-Newton and NuSTAR. Moreover, the two combined HEX-P instruments yield high sensitivity to pulsations in both soft and hard sources, something difficult to achieve with current instruments without complicated coordination. Finally, the LET’s short frame time allows for sensitive searches of pulsations at frequencies well above XMM-Newton/EPIC-pn’s limit.

FIGURE 4. Results of a blind search of pulsations described in Section 3.2.3. Lines connect candidates for the same pulsar in different instruments. Some patterns appear: 1) XMM-Newton/EPIC-pn and HEX-P/LET have similar performance at low frequencies (the balance depends on pn’s slightly larger effective area and LET’s lower contamination), while HET is always better than NuSTAR; 2) EPIC-pn loses sensitivity above $\sim 7$ Hz due to its frame time, and LET above ∼200 Hz.

To measure the expected detection significance of pulsars injected in the SIXTE simulation, we first evaluated the H-test at the injected frequency and the integrated pulse profile over the full simulation, and found the correct number of harmonics M to describe the pulsed profile. We then estimated the variation in the $Z_M^2$ statistic for different exposures and source distances in the following way:

1. We created a grid of exposure times (E) and distances (D) for the source. For each value of E and D, we rescaled all sky counts (background and source) with a E/D² law, while the instrumental background was rescaled by exposure alone.

2. For each of the rescaled count estimates, we constructed a model pulse profile, as follows: 1) the theoretical pulse profile, including its intrinsic pulsed fraction, was normalized to an integrated value of 1 and multiplied by the number of source counts; 2) we added a background with the expected counts per bin from the sky and the instrument.

3. Finally, we created 100 realizations of the pulse profile by simulating Poisson-distributed random numbers for each bin of the pulse profile, and we calculated the $Z_M^2$ statistic for each realization. We then calculated the false-alarm probability for the average value of the $Z_M^2$ statistic [using Buccheri et al. (1983), as adapted to folded pulsed profiles by Bachetti et al. (2021)] over all realizations, and expressed it in terms of Gaussian sigma values for ease of interpretation. This was then used to form the plots on the right hand side in Figure 5.

Figure 5 shows the detectability of two simulated bright pulsars (not even at a ULX level of brightness), near the nucleus of a galaxy modeled on NGC 253, with different instruments; these provide good examples of what happens when one looks for pulsations from sources in a crowded field and were selected to have a softer and harder spectrum respectively. The detectability of specific pulsars depends on the spectral shape, the fraction of pulsed flux, and the sky and instrumental background. Therefore, the same intrinsic pulsed flux ratio (PFR)⁶ can result in very different total PFR (and hence, detectability) in different instruments, depending on sensitivity and angular resolution. In these two cases (and the same consideration is supported by the blind search as well, Figure 4), one can see why the instruments onboard HEX-P perform equivalently or better than current focussing instruments like XMM-Newton and NuSTAR, especially if the source is located in crowded fields or even in peripheral regions of the host galaxy. The angular resolution of HEX-P allows the emission from a large number of contaminating sources to be filtered out, improving the PFR and the pulsation sensitivity, even where the effective area is slightly lower (e.g., XMM-Newton/EPIC-pn compared to LET).

FIGURE 5. The maps on the left show a 1-Msec LET simulation of the central region of NGC 253 using the same setup as Figure 1. Each map is centered on one example source in this region, and the two insets show the spectrum and the pulse profile of the source. Circles indicate the HPD of each instrument. The two sources are only resolved by the LET. The shaded regions on the right of each map show (top) the limiting distance at which we can detect the pulsation (at 5 σ) at the given PFR with different exposures, and (bottom) the limiting intrinsic PFR detectable at the distance of NGC 253 (3.5 Mpc, marked in the top plot as a blue dashed line). Color coding for all line plots indicates different instruments (the band pass in the first inset). NuSTAR and HEX-P/HET include the flux from two modules. Percentages in the legend indicate the source to background (all components) flux ratio in the region.

3.3 Spectral simulations

As discussed in Section 1, the X-ray spectral continuum of ULXs can be complex, characterized by up to three main components: a low-temperature ($<1$ keV) thermal component, a hotter 1–3 keV component, which appears broader than the cool one, and a hard component above 10 keV (e.g., Walton et al., 2018b). The latter can appear as a hard excess in the NuSTAR band, and may be related to Compton up-scattering of the disk emission, which is not included in the simple accretion disk models typically used [c.f. post-processing of general relativistic magnetohydrodynamical simulations of super-Eddington discs, (Mills et al., 2023)]. In some cases it is likely associated with the accretion column of a NS. The cool thermal component is often modelled with a simple disk-blackbody component (e.g., diskbb in Xspec) approximating the emission from the outer disk and/or wind. The hotter thermal component has a broader spectrum and can be equally described by a slim disk (e.g., diskpbb in Xspec) with a temperature-radius law T ∝ R^−p (with p taking values 0.5–0.75, e.g., Walton et al., 2020) or a thermal Comptonisation component (e.g., nthcomp in Xspec, Middleton et al., 2015a). Finally, when identified with a pulsation, the emission from the accretion column can be empirically modelled with a cutoff-powerlaw (cutoffpl) with slope Γ = 0.59 and high-energy cutoff E_cut = 7.9 keV, which is based on the average parameters found for known ULX pulsars via phase-resolved spectroscopy (e.g., Walton et al., 2018a).

The study of the relationship between the luminosities and temperatures of the thermal components measured in different epochs may be useful in order to obtain valuable information on the evolution and structure of the inner accretion disk. The luminosity-temperature (L − T) trends often disagree with the theoretical predictions of sub-Eddington thin discs (L ∝ T⁴) and even show negative slopes consistent with the photosphere of a super-Eddington disk/wind (e.g., Robba et al., 2021). In this regard, HEX-P will provide an important improvement with respect to current missions thanks to its high effective area and broad energy coverage. In Walton et al. (2020) (in prep.), for instance, simulations of Holmberg IX X-1 have been performed showcasing how HEX-P will enable accurate estimates of the L − T trends. In particular, it will be possible to distinguish the L − T tracks observed at high and low source fluxes [also seen in NGC 1313 X-1 (Walton et al., 2020)], which can indicate dramatic changes in the inner accretion flow.

Once a satisfactory description of the continuum is achieved, spectral residuals can be located [almost ubiquitously: Middleton et al. (2015b)] and, eventually, searches for narrow lines can be performed. This is typically performed through scans that involve a moving Gaussian through a grid of centroid energies (e.g., Pinto et al., 2016). Monte Carlo simulations must be used in order to probe the statistical significance of any putative lines (e.g., Kosec et al., 2021). More recently, the exploration of a large parameter space has been employed through grids of plasma models, which account for multiple emission or absorption lines simultaneously. This boosted the detectability of weak lines such as those from ultrafast outflows (e.g., Kosec et al., 2018; Pinto et al., 2021). The emission lines are primarily produced by Ne K and Fe L ions around 1 keV and the absorption lines from blueshifted O VIII and Ne K/Fe L around 0.8 and 1.2 keV, respectively. Evidence of Fe K around 8–9 keV has also been found, showing the hottest, fastest and most powerful wind component (e.g., Walton et al., 2016b; Brightman et al., 2022). The emission and absorption lines can be well described with models of plasma in photoionization equilibrium, which is expected in the case of winds irradiated by the ULX continuum. Most work in this regard has made use of the Spex code (Kaastra et al., 1996), in particular through the pion component, although some work also made use of Xstar (Kallman and Bautista, 2001; Mendoza et al., 2021) coupled with Xspec (Arnaud, 1996). In the latter case, photoionization model grids are computed with Xstar and then loaded into Xspec, whilst in the former, Spex instantaneously fits the continuum, calculates the ionization balance and the model, and fits it to both lines and continuum. An example of HEX-P’s performance with regards to continuum modelling and line detection is shown in Section 3.3.2.

3.3.1 Simulation setup

A typical workflow for X-ray spectral simulations involves: 1) defining a theoretical model; 2) defining a count rate $\mathcal{C}$ and observation length $\mathcal{L}$; 3) specifying a set of instrument responses (RMF and ARF); 4) specifying an appropriate background file or a model providing an estimate of the instrumental and sky background contaminating the source region; 5) using a tool that generates a list of $\mathcal{C}\mathcal{L}$⁷ random X-ray energies following a statistical distribution obtained by the convolution of the spectral model with the instrumental responses.

In our case, the last step in the above list is executed with Xspec’s fakeit command, using the most recent response files available for each mission (see Sections 3.3.2, 3.3.3), following the procedure from the Xspec manual⁸ or the equivalent procedure for PyXspec⁹. In order to simulate emission or absorption lines from plasmas in different types of equilibria (collisional or photoionisation) we also used the Spex code (Kaastra et al., 1996), which provides up-to-date atomic data, plasma models and a fitting package¹⁰.

3.3.2 Wind features

In order to demonstrate the power of HEX-P in studying winds from ULXs by, e.g., achieving a highly significant $(\sim 5\sigma )$ Fe K detection, we performed a simulation of the archetypal ULX NGC 1313 X-1, adopting the model described in Pinto et al. (2020b). The continuum consists of cool (0.2 keV) and hot (2 keV) thermal components and a cut-off powerlaw, the latter dominating above 10 keV with a slope of Γ = 0.59 and high-energy cutoff E_cut = 7.9 keV as mentioned above (Walton et al., 2020). Starting from previous work on NGC 1313 X-1, we simulated a deep HEX-P observation adopting the best-fit multiphase wind model of photoionised emission (at rest) and absorption (blueshifted by 0.2c) plasmas [see Pinto et al. (2020b)]. This consists of one PION emission component and two XABS absorption components in Spex. The results are shown in Figure 6. The LET and HET simulated spectra in the top panel show a very strong absorption feature in the soft band above 1 keV from the cooler phase of the wind, and a narrow absorption line just above 8 keV from the hot phase of the wind [mainly Fe XXV; c.f., (Walton et al., 2016b)]. We also highlight the contribution of the wind by removing the emission and absorption lines from the model and refit the simulated spectra with a continuum-only model [double thermal component plus cutoff powerlaw; see Walton et al. (2020)]. Each of the three emission and absorption line components is detected¹¹ at ≳5σ assuming an exposure of 500 ks, whilst for the weak Fe K lines confidence levels around 3 σ are found with 300 ks observations, which are significantly shorter than the 700 ks total XMM-Newton and NuSTAR observations needed. Moreover, any observation that can detect the Fe K absorption will easily detect the lower energy ($\sim 1$ keV) absorption as well since the latter typically requires exposures of a few tens of ks for bright ULXs. There are about 10 ULXs with fluence comparable to NGC 1313 X-1 in the Fe K band and a few dozen within a factor of a few (Walton et al., 2022). We therefore be confident that HEX-P will enable the first meaningful search for Fe K lines in a statistical sample of ULXs, which will place new constraints on the wind geometry, its outflow and kinetic energy, and improve our understanding of the overall accretion mechanism.

FIGURE 6. 500ks HEX-P simulation of ULX NGC 1313 X-1 adopting the best-fit multiphase wind model of photoionised emitting and absorbing plasma (Pinto et al., 2020b). The middle panel shows the residuals assuming a continuum-only model (double blackbody plus cutoff powerlaw). The bottom panel shows the residuals to the full wind model. The rest-frame energies of the dominant transitions are labelled. Each component producing emission/absorption lines is detected at ≳5σ.

3.3.3 Cyclotron lines

Three (potential) CRSF features have been discussed for ULXs/related sources in the literature: an unusually low-energy ∼4.5 keV feature in M51 ULX-8 that could be a proton CRSF [seen by Chandra; (Brightman et al., 2018)], an extremely high-energy feature at ∼150 keV in the Galactic, Swift J0243+6124 [seen by HXMT; (Kong et al., 2022)], and a potential CRSF at ∼13 keV in the ULX pulsar NGC300 ULX1 [seen in a coordinated XMM-Newton + NuSTAR observation; (Walton et al., 2018a)]. However, the detection of the latter is noted to be somewhat tentative, as the data can also be explained by a more complex continuum model with no CRSF, an issue also discussed by Koliopanos et al. (2019) who argued the more complex continuum model may actually be preferred by the data. Part of the uncertainty is the poorly understood continuum and the limited coverage above ∼30 keV offered by NuSTAR, which makes it an excellent case-study for HEX-P.

In order to estimate the exposure time needed to unambiguously detect the putative line in NGC 300 ULX-1, we have simulated spectral models with and without the feature, using the models reported by Koliopanos et al. (2019) (see their Table 1 where they use a combination of two multi-color disk blackbodies and a powerlaw for the high energy emission, referred to as their MCAE model). We focus our simulations on the MCAE model with and without the line, as this gave the highest (Bayesian) evidence among a set of competing models (see their Table 2), but below we also present an alternative case.

We estimated the CRSF detection significance α, using the method from Protassov et al. (2002). To this end, we simulated a single HEX-P spectrum (both LET and HET) using the MCAE model with a CRSF line [as reported in Koliopanos et al. (2019)] for exposure times of 50, 100, 150 and 200 ks. We fitted this spectrum with the null hypothesis model of the MCAE model without a line and then with the MCAE model with the CRSF, recording the single value of $\mathrm{\Delta }{\chi }_{\text{CRSF}}^2$ fit-improvement for each exposure time. We then simulated 1,000 spectra from the best-fit null hypothesis model and fitted these with both models, also recording the fit-improvement in each instance (obtaining the distribution of $\mathrm{\Delta }{\chi }_{\text{null}}^2$). We compared this fit-improvement distribution $\mathrm{\Delta }{\chi }_{\text{null}}^2$ with the $\mathrm{\Delta }{\chi }_{\text{CRSF}}^2$ value obtained from the HEX-P spectrum that included the CRSF, thereby obtaining the probability of observing a Δχ² fit-improvement as high as that observed in the original HEX-P spectrum (i.e., the line detection significance α). We repeated the whole process 10 times per exposure to estimate the variations in $\mathrm{\Delta }{\chi }_{\text{CRSF}}^2$ and α due to Poisson fluctuations only.

The results of our simulations are shown in Figure 7. We see that for a 50 ks exposure, including the line represents a $\mathrm{\Delta }{\chi }_{\text{CRSF}}^2\sim 19\pm 6$ the 1σ uncertainty comes from the spread of the 10 simulations improvement over the null continuum, which translates to a ∼97% (2.2σ) detection significance, although the high-energy power-law will be marginally detected. Such a $\mathrm{\Delta }{\chi }_{\text{CRSF}}^2$ is already higher than, or comparable to, the value obtained by Koliopanos et al. (2019), who reported a $\mathrm{\Delta }{\chi }_{\text{CRSF}}^2=12$ improvement based on the existing ∼140 ks EPIC-pn and 180 ks NuSTAR exposure (see Figure 7 for details). For a 100 ks exposure time with HEX-P, the line is detected at the $\sim 99\%$ $(\sim 2.6\sigma )$ confidence level. Here the line provides a $\mathrm{\Delta }{\chi }_{\text{CRSF}}^2$ = 29 ± 10 fit-improvement. Finally for exposure times greater than 150 ks, we found the line will be detected comfortably above the 99.73% (3σ) level with $\mathrm{\Delta }{\chi }_{\text{CRSF}}^2\sim 50\pm 14$ and 58 ± 14 for 150 and 200 ks respectively, with the high-energy power-law significantly detected. This final spectrum as seen by HEX-P is shown in Figure 8.

FIGURE 7. HEX-P detection levels and constraints on the putative NGC 300 ULX-1 CRSF. (Left) Line detection significance (top) and Δχ² fit-improvement (bottom) as a function of exposure time. The error bars represent the spread over 10 realizations (see text for details). The Δχ² fit-improvement obtained by Koliopanos et al. (2019) using a combined XMM-Newton and (∼140 ks for EPIC-pn) NuSTAR (∼180 ks) observation is shown as a dashed line for comparison. (Right) Constraints on the NGC 300 ULX-1 CRSF centroid energy (top) width (middle) and depth (bottom) as a function of exposure time for HEX-P (orange) and XMM-Newton + NuSTAR (blue, values have been shifted slightly horizontally for clarity). The mean best-fit value and mean (90% level) uncertainty are shown after averaging over 25 spectra per exposure time. The stars indicate the constraints obtained by Koliopanos et al. (2019) for comparison.

FIGURE 8. (Top) CRSF line in NGC 300 ULX-1 as observed by HEX-P in a 200 ks observation. (Bottom) CRSF line in M51 ULX-8 as observed by HEX-P in a 100 ks observation.

We also obtain simulated constraints on the line energy, width and depth, and compare these to those reported by Koliopanos et al. (2019) in Figure 7. We also simulated XMM-Newton-NuSTAR observations including all EPIC and both FMPA/B cameras, using responses from a similar epoch to the NGC 300 ULX-1 observations. We note that our latter simulations ignore calibration inaccuracies between XMM-Newton and NuSTAR detectors, which will only degrade the results presented here, particularly as compared to HEX-P. We can see from Figure 7 that HEX-P will allow us to obtain tighter constraints on the line properties compared to dedicated XMM-Newton + NuSTAR observations (about a factor two or even more). As stated above, one of the main uncertainties in the identification of the CRSF in NGC 300 ULX-1 was the nature of the underlying continuum. Walton et al. (2018a) proposed to tackle the problem by performing pulse-phase resolved spectroscopy. In particular, Walton et al. (2018a) isolated the pulsed spectrum by extracting data from the brightest and faintest quarter-phases of the pulse cycle of NGC 300 ULX-1 (and obtained the difference spectrum, i.e., “pulse-on”—“pulse-off”). By doing so, the uncertainty on the continuum is reduced as the ‘constant’ or non-pulsing component is eliminated. We thus performed additional simulations to examine the ability of HEX-P to distinguish between the two different continuum models presented by Walton et al. (2018a).

We simulated separate (LET and HET) spectra for the high (pulse on) and low (pulse off) phases of the pulse cycle of NGC 300 ULX-1 (each spanning 0.25 in phase) using the models of Walton et al. (2018a), i.e., a combination of two non-pulsing blackbodies and a high-energy gabs⊗cutoffpl ascribed to the accretion column. We subtracted the two and then fitted the pulsed component with the two competing models presented by Walton et al. (2018a) (a model containing the CRSF line, gabs ⊗ cutoffpl, and a more complex featureless continuum, cutoffpl ⊗ simpl) and retrieved the difference in Bayesian Information Criterion [ΔBIC; (Schwarz, 1978)] between the two models for 100 simulations per exposure time. From Figure 9 we can see that HEX-P would be about a factor 1.5 more effective (notably above ∼50 ks per phase bin, corresponding to a total 200 ks) than XMM-Newton & NuSTAR in distinguishing different continuum models, therefore reducing uncertainties related to the presence of CRSF lines in the ULX spectra.

FIGURE 9. ΔBIC between the two competing models presented in Walton et al. (2018a) (a simpler continuum including the CRSF vs. a more complex featureless continuum) fitted to the pulsed component of NGC 3000 ULX1 for HEX-P (top) and combined XMM-Newton + NuSTAR observations (bottom) for different exposure times. 100 realizations are shown for each exposure time, as indicated by the colors. For comparison, we also show the value retrieved by Walton et al. (2018a) using 140 of EPIC-pn, 190 ks of combined EPIC-MOS, and 180 ks of combined NuSTAR observations (dashed line). We can see that HEX-P would be about a factor 1.5 more effective (notably above ∼30 ks) in distinguishing these two different continuum models. Note that the exposure times are given for a quarter of the full pulse cycle and therefore total exposure times would be four times the values shown.

The same procedure as above was followed to explore the ability of HEX-P to detect the 4.5 keV CRSF originally detected in the Chandra spectrum of M51 ULX-8 by Brightman et al. (2018). We used the models reported by Middleton et al. (2019a), who carried out a more detailed analysis of the continuum (modelled with a soft diskbb and a cutoffpl for the high energy emission). We found that a 50 ks HEX-P observation detects the line at the ∼97% (2.2σ) level while a 100 ks exposure detects the line in excess of 99.73% (3σ) (to be compared with the 99.98% or 3.8σ detection by Chandra in 180 ks). Figure 8 shows the spectrum as observed by HEX-P.

Another unsolved problem HEX-P will tackle is determining the population of particles producing the CSRF lines. In the case of electrons, further associated lines at the harmonic frequencies are expected, whereas protons are expected to produce lines only at the fundamental frequency [see Staubert et al. (2019) and references therein]. Brightman et al. (2018) were unable to detect any harmonic feature at 9 keV, likely due to the low effective area of Chandra in this band. If we assume an additional line with the same properties (width and depth) as the fundamental but at twice the energy, the significance of the CRSF detection (now combined) would obviously increase (for 100 ks, an additional ∼Δχ² = 26 compared to the case where only the fundamental is present in the spectrum). We found that an exposure of about 800 ks using HEX-P would allow us to detect the harmonic at the 99% (2.6σ) level alone. At this significance, combining both features will allow us to unambiguously identifying whether the CRSF is due to electrons or protons and thus accurately determine the magnetic field strength close to the NS.

4 Conclusion

As extragalactic sources, every ULX study–population, spectroscopy, temporal analysis–is affected by the sensitivity, angular resolution and instrumental background of X-ray telescopes. All major discoveries in the field of ULXs have been driven by instruments combining focusing capabilities and effective area over the $\sim 0.3-30$ keV band, such as the combination of XMM-Newton and NuSTAR. However, many galaxies containing ULXs are still unresolved by these two instruments, e.g., M51 and NGC 253. In addition, the effectiveness of multi-instrument analysis is lowered by the difficulties of retrieving simultaneous observations.

HEX-P’s low- and high-energy telescopes ensure broadband simultaneous coverage of the X-ray spectrum and excellent angular resolution, significantly improved over both XMM-Newton (for the LET) and NuSTAR (for the HETs). Thanks to the lower background ensured by its small point-spread function, HEX-P facilitates the identification of ULXs through their distinct curved spectra even at large distances.

The design of the HEX-P mission would make it the ideal facility to shed light on a large number of questions regarding ULXs and super-Eddington accretion in general. In particular, broad-band and good energy resolution spectra would lead to a deep understanding of the balance between the mass inflow and outflow. HEX-P will also enrich ULX population studies, for example, by identifying the presence of NSs in ULXs through the detection of pulsations and CRSFs with significantly higher fidelity than existing missions. As shown in Section 3.2, HEX-P surpasses the sensitivity of NuSTAR and XMM-Newton in detecting pulsations from sources in crowded fields, especially ULXs. As a result, HEX-P extends the range for detecting pulsars to distances 2–5 times farther than currently achievable, thereby probing a much larger (8–125 times) volume. This will permit access to a statistically significant population of bright extragalactic X-ray pulsars for the very first time. Further demonstration of HEX-P’s capabilities in pulsar searches in crowded environments can be found in Mori et al. (Galactic Center) and Alford et al. (2023) (Magnetars, CCOs) in this volume.

In addition to the above, the L1 orbit of HEX-P enables a continuous observation of a large fraction of the sky, without Earth occultation (an issue for NuSTAR, which is located in a low Earth orbit) or increased particle background from the Van Allen belts (an issue for XMM-Newton and Chandra, with their highly elliptical orbits). This allows a significantly shorter “clock time” (actual time spent observing, as opposed to clean source exposure), while the broad bandpass of HEX-P eliminates the issues due to differing observation windows when coordinating satellites in different orbits (Middleton et al., 2017). This makes HEX-P more efficient, enabling observations of more sources in a given amount of time, with much longer exposure times for observations constrained by clock time (e.g., coordinated observations, orbital measurements).

In summary, the HEX-P mission concept holds tremendous promise for advancing our understanding of ULXs and the physics behind super-Eddington accretion. Its broad energy coverage, improved sensitivity, and exceptional spectral and timing capabilities establish it as a powerful tool for investigating ULXs and the fundamental physics at play.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

MaB: Conceptualization, Formal Analysis, Investigation, Software, Supervision, Visualization, Writing–original draft, Writing–review and editing. MM: Conceptualization, Supervision, Writing–original draft, Writing–review and editing. CP: Formal Analysis, Investigation, Methodology, Visualization, Writing–original draft. AG: Formal Analysis, Investigation, Methodology, Visualization, Writing–original draft. DW: Conceptualization, Supervision, Writing–review and editing. MuB: Writing–review and editing, Conceptualization. BL: Writing–review and editing, Methodology, Software. TR: Writing–review and editing, Conceptualization. GV: Writing–review and editing, Conceptualization, Investigation, Methodology. JA: Writing–review and editing. RA: Writing–review and editing. EA: Writing–review and editing. LD: Writing–review and editing. HEa: Writing–review and editing. HEl: Writing–review and editing. JG: Writing–review and editing, Funding acquisition, Supervision. GL: Writing–review and editing. AJ: Writing–review and editing. KMa: Writing–review and editing, Funding acquisition, Supervision. CM: Writing–review and editing. SM: Writing–review and editing. KMo: Writing–review and editing, Software. FP: Writing–review and editing. KO: Writing–review and editing. MP: Writing–review and editing. DS: Writing–review and editing, Funding acquisition, Supervision. GY: Writing–review and editing. AW: Writing–review and editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. MaB was funded in part by PRIN TEC INAF 2019 “SpecTemPolar!—Timing analysis in the era of high-throughput photon detectors.” TR acknowledges funding from STFC as part of the consolidated grants ST/T000244/1 and ST/X001075/1. GV acknowledges support by Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “3rd Call for H.F.R.I. Research Projects to support Postdoctoral Researchers” through the project ASTRAPE (Project ID 7802). The work of DS was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

¹See Ludlam et al. (2023), this volume, for CRSF studies on bright Galactic sources with HEX-P.

²They can be downloaded at https://www.sternwarte.uni-erlangen.de/sixte/instruments/

³The galaxy does not have a bright AGN in the center. Such a source could, in some cases, dominate the emission close to the nucleus, and would be a source of confusion for the detection of ULXs for all instruments, plausibly more so for XMM-Newton and NuSTAR, that have a broader PSF.

⁴See, e.g., the XMM-Newton science analysis software documentation, https://xmm-tools.cosmos.esa.int/external/xmm_user_support/documentation/sas_usg/USG/

⁵Here, we are neglecting data loss due to occultation or flaring during the orbit. Taking them into account would advantage HEX-P further, because the same exposure would be obtained in shorter observations.

⁶Note that with PFR we are comparing the pulsed flux, the integral of the AC component of the pulse, to the total flux. The more commonly used pulsed fraction (PF), instead, compares the amplitude of the pulsation to the mean or the maximum of the pulse. Our definition is more independent of pulse shape.

⁷Or better, a random number following a Poisson distribution around that value.

⁸https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/node41.html

⁹https://heasarc.gsfc.nasa.gov/xanadu/xspec/python/html/extended.html

¹⁰https://spex-xray.github.io/spex-help/index.html

¹¹As determined through Δχ2, using the affected spectral bins as the number of free trials.

References