© The Authors 2025

1. Introduction

The number of ultraluminous X-ray sources (ULXs) comprises about two thousand objects (e.g. Walton et al. 2022) characterised by high assumed isotropic X-ray luminosities well in excess of the Eddington limit for a 10 M_⊙ black hole (BH; for recent reviews, see e.g. Kaaret et al. 2017; King et al. 2023; Pinto & Walton 2023). Although some of the ULXs are still believed to host sub-Eddington accreting intermediate mass BHs (10² − 10⁵ M_⊙; Colbert & Mushotzky 1999), most of the population is likely composed by super-Eddington accreting stellar mass compact objects, namely, BHs (M_BH = 5 − 80 M_⊙) or neutron stars (NSs). The number of confirmed NSs in the ULX population has grown significantly in the past decade, with the discovery of at least six (extragalactic) pulsating ULXs (PULXs), which is the smoking gun for an NS nature of a compact object in such sources (Bachetti et al. 2014; Israel et al. 2017a; Carpano et al. 2018; Israel et al. 2017b; Fürst et al. 2017; Sathyaprakash et al. 2019; Rodríguez Castillo et al. 2020). Luminosities in PULXs can be very extreme, in some cases up to 1000 times their Eddington limit (Israel et al. 2017a). PULXs reside in binary systems with young and massive companion stars where orbital periods range between 1 and 100 days and the projected semi-major axis is of a few to tens of light seconds (e.g. Bachetti et al. 2022; Belfiore et al. 2024). Pulsations in PULXs are intermittent and characterised by a small pulsed fraction, which makes increasing the number of NSs in ULXs challenging. However, it cannot be excluded that the NS population in ULXs is higher than currently known. On the basis of the PULX properties, indirect pieces of evidence (e.g. Middleton & King 2017; King & Lasota 2020) can be adopted to unveil candidate NSs in other ULXs. In particular, this evidence can be based on a switch between accretion and quiescent epochs during the ULX lifetime, which is believed to be due to the ignition of the propeller (e.g. Tsygankov et al. 2016; Earnshaw et al. 2019; Song et al. 2020; Chashkina et al. 2017; Middleton et al. 2023); a hard spectrum (e.g. Pintore et al. 2017; Koliopanos et al. 2017; Walton et al. 2018, 2020; Gúrpide et al. 2021), likely associated with the emission from an accretion column at the NS polar caps; or the detection of cyclotron absorption lines as in the case for M51 ULX-8; Brightman et al. 2018; Middleton et al. 2019).

In order to achieve luminosities as high as to 10⁴¹ erg s⁻¹, a stellar-mass compact object needs to accrete matter well above the Eddington limit (L_Edd = 1.4 × 10³⁸ (M/M_⊙) erg s⁻¹). At accretion rates Ṁ ≳ Ṁ_Edd = L_Edd/(ηc²), where η is the efficiency and c the speed of light, radiation pressure is expected to thicken the accretion disc, which then exhibits thermal X-ray spectra peaking at a few keV with a high-energy cut-off around 10 keV. This has been confirmed by X-ray broadband ULX spectra obtained with XMM-Newton and NuSTAR (e.g. Bachetti et al. 2013; Walton et al. 2014, 2020). Moreover, above the critical rate Ṁ_c = 9/4 Ṁ_Edd (Poutanen et al. 2007), powerful outflows of ionised plasmas are expected to be driven by radiation pressure, as shown by magneto-hydrodynamical simulations (e.g. Ohsuga & Mineshige 2011). These winds imprint weak lines in the high-resolution X-ray spectra of ULXs (e.g. Pinto et al. 2016; Kosec et al. 2021, for a review see also Pinto & Kosec 2023), primarily in the soft X-ray band (0.3–2 keV). The discovery of winds in two pulsating ULXs, NGC 300 ULX-1 (Kosec et al. 2018a) and the Galactic Swift J0243.6+6124 (van den Eijnden et al. 2019), implies that the dipole component of the magnetic field is ≲10¹³ G; otherwise, the latter would truncate the disc beyond the spherisation radius (R_sph = 27/4 Ṁ R_G) and prevent radiation pressure from launching the observed, relativistic (0.22c), wind. Although inclination effects, ionisation of the accretion flow, or beaming can be non-negligible, the discovery of a wind with a certain escape velocity might therefore place constraints on the field strength (for the dipolar component) and the nature of the accretor.

The ULX NGC 4559 X7 (hereafter X7), in the galaxy NGC 4559 (distance 7.5 Mpc; Tully et al. 2016), is a peculiar object with a marked flaring activity during the epochs of its highest observed flux states (Pintore et al. 2021). Such flares were shown to be driven by flux changes at energies higher than 1 keV, with a clear harder-when-brighter behaviour (Pintore et al. 2021), which is seen in several ULXs (e.g. Gúrpide et al. 2021). High-resolution spectroscopy data from XMM-Newton/Reflection Grating Spectrometer (RGS) data have revealed a number of narrow absorption and emission lines that are compatible with the existence of an outflow in X7. In addition, variations in the hard band (2 − 10 keV) were seen to lag behind the soft band (0.3 − 2 keV) with a delay of hundreds of seconds (e.g. De Marco et al. 2013; Pinto et al. 2017; Kara et al. 2020), suggesting that we may be seeing accretion rate fluctuations propagating inwards through the accretion flow. However, the nature of the compact object powering X7 still remains puzzling.

In this work, we extend the analysis shown in Pintore et al. (2021), and we further investigate the short- and long-term spectral and temporal properties of X7. Overall, our aim is to identify if the system hides a super-Eddington accreting BH or an NS.

2. Data reduction

2.1. XMM-Newton

XMM-Newton observed NGC 4559 on 27 May 2003 for ∼42 ks (hereafter XMM0; ObsID: 0152170501, PI: M. Cropper). Then our group obtained a new ∼75 ks-long observation (hereafter XMM1; ObsID: 0842340201, PI: F. Pintore), which was taken on 16 June 2019, and a third one on 13 June 2022 of ∼130 ks (hereafter XMM2; ObsID: 0883960201, PI: C. Pinto). XMM0 and XMM1 were already reported in Pintore et al. (2021). XMM1 and XMM2 caught the source during its highest fluxes (see the next sections).

We reduced the observations with SAS v20.0.0, using data from both EPIC-pn and EPIC-MOS (1 and 2) cameras: events were selected with PATTERN ≤ 4 (i.e. single and double-pixel events) for the pn, and PATTERN ≤ 12 (i.e. single- and quadruple-pixel events) for the MOS. All observations are affected by high background time intervals, which were removed from the analysis. This resulted in a net exposure time of ∼35 ks, ∼65 ks and ∼95 ks for XMM0, XMM1 and XMM2, respectively. We extracted source and background events from circular regions with radii of 30″and 60″, respectively. The photon times of arrival (ToAs) were converted to the barycentre of the Solar System with the task BARYCEN, using the best Chandra coordinates of the target (RA = 12h 35m 51.71s, Dec = +27d 56m 04.1s; Swartz et al. 2011) and ephemeris DE405¹.

EPIC spectra were grouped with the FTOOLS GRPPHA in order to accumulate at least 25 counts per energy bin. The RMFGEN and ARFGEN tools were used to generate response matrices and ancillary files, respectively. For all the EPIC light curves, we applied the SAS task EPICLCCORR to take into account vignetting, point-spread-function, quantum efficiency and bad pixels.

High-resolution spectra require long exposures to estimate the background accurately and to fill the several thousand energy channels. We focussed on the data from the RGS of both the XMM1 and XMM2 observations only, because the short exposure of XMM0 was not enough for RGS analysis of its data. The RGS data were reduced with the RGSPROC task, which produces calibrated event files, spectra and response matrices. Following the standard procedures, we filtered the RGS data for solar flares using the background light curve from the RGS CCD number 9 (corresponding to ≲7.5 Å or ≳1.7 keV). The background light curves were binned at a 100 s time resolution and all the time bins with a count rate above 0.2 cts s⁻¹ were rejected,

using the same good time intervals for the RGS 1 and 2 detectors. We extracted both the first- and second-order RGS spectra in a cross-dispersion region of 1 arcmin width, centred on the same source coordinates used for extracting the EPIC spectra. The background spectra were extracted from regions outside the 98% of the source point-spread-function. Due to the possible presence of further fainter sources in the background regions, we also extracted the background spectrum produced using the template background model from blank field data (obtained from the rgsproc pipeline with the flag ‘withbackgroundmodel=yes’). The two background spectra were consistent. In order to obtain the highest signal-to-noise ratio, we stack the first- and second-order spectra from the two observations for a total of 176 ks for each spectral order. Such a combined RGS spectrum provides about 9400 net counts, which broadly agrees with the ≳10 000 counts typically needed to significantly detect lines in ULXs (Pinto et al. 2016; Kosec et al. 2021). To search for spectral features, we also need an accurate knowledge of the continuum shape, especially for photoionisation modelling. We therefore stacked the MOS 1, MOS 2, and pn spectra of XMM1 and XMM2 to cover the 2–10 keV band. In the end, we have 5 combined spectra for the XMM1 and XMM2 observations (RGS first and second order, MOS1, MOS2, and pn).

2.2. Swift

The X-Ray Telescope (XRT) instrument onboard the Neil Gehrels Swift Observatory (hereafter Swift/XRT) observed the galaxy NGC 4559 for 164 times between January 2006 and July 2024 (Log in Table B.1). We note that the vast majority of the observations were taken in photon counting (PC) mode with an average exposure time of about 2 ks each, generally taken with a weekly or bi-weekly cadence. We excluded all the observations in other modes (e.g. windowed timing) in order to avoid contamination from other sources in the field of view.

We reduced all observations with the XRTPIPELINE tool of the FTOOLS software, and we extracted source and background events from circular regions of radii 40″and 60″(centred at the source position and in a source-free region of the CCD), respectively. The long-term light curve of the Swift/XRT observations is shown in Fig. 1. For a comparison with the XMM-Newton, NuSTAR and Chandra data, we converted Swift/XRT count rates in luminosities. Due to the low statistics of the single Swift/XRT observations, we performed the count-rate to flux conversion by firstly fitting the stacked spectrum² of all Swift/XRT observations with an absorbed power-law spectral models. We determined the best-fit by fixing the column density to 1.5 × 10²¹ cm⁻², which is the value found from the high quality XMM-Newton spectra (see Sect. 3.3.1), and let the power-law photon index Γ and normalisation free to vary. The best fit of the stacked spectrum gives Γ ∼ 2.2 and 0.3–10 keV flux of 1.1 × 10¹² erg cm⁻² s⁻¹. We adopted this model, freezing the photon index and letting free to vary only the normalisation, to estimate the observed 0.3–10 keV source luminosity in each Swift/XRT observation.

Fig. 1. Light curve of X7 in the 0.3–10 keV energy band when using all available X-ray observations. The XMM-Newton, NuSTAR, Chandra, and Swift/XRT luminosities are not corrected for absorption and are taken as the average in each observation (assuming a distance of 7.5 Mpc). Chandra data are taken from Pintore et al. (2021). The inset shows the light curve after 2018 with the super-position of the profile of the tentative ∼240 d super-orbital periodicity (magenta).

3. Results

3.1. Long-term analysis

As reported in Pintore et al. (2021), a tentative super-orbital periodicity of ∼200 d was proposed to describe the (still limited at that time) Swift/XRT long-term light curve of X7. Thanks to the regular Swift/XRT monitoring performed between 2020 and 2024, the baseline to search for periodicities is significantly increased.

We ran a Lomb-Scargle (LS) analysis over all observations between 2020–2024 and investigated the energy bands 0.3–10 keV, 0.3–1.0 keV and 1.0–10 keV, with a bin time of 2 days. We used the LombScargle class of the python Astropy library, where we set the method to initially estimate the false alarm probability³ of any signal to Baluev (Baluev 2008). The LS, applied separately to the three energy bands, revealed an indication for a periodicity in the range 500–600 d with a 3σ significance in the total and the 0.3–1.0 keV bands. We note that the peak is very broad, pointing towards large uncertainties. In addition, such a periodicity does not cover more than 2–3 cycles of the whole light curve, suggesting that it should be only treated as a possible spurious signal due to the red noise.

On the other hand, we note that excluding the observations after 21 January 2023 (MJD 59965), i.e. during which the source showed a prolonged and stable low flux, the LS provides again the periodicity of ∼500 − 600 d but also signals at 239 d (significance > 3σ), 238 d (significance > 3σ), and 242 d (slightly less than 3σ) in the 0.3–10 keV, 1–10 keV and 0.3–1.0 keV energy bands, respectively (Fig. 2-top).

Fig. 2. Top: Lomb-Scargle periodograms of the Swift/XRT observations in the time intervals between 2020 and January 2023 in the 0.3–10 keV, 0.3–1 keV, and 1–10 keV energy range. Dashed lines indicate signals with > 3σ significance. Bottom: Full 1–10 keV Swift/XRT light curve folded with a period of 240 d.

We assess the uncertainty on the ∼240 d periodicity by performing a Monte Carlo simulation. For each point of the real Swift/XRT light curve, we sampled a random value following a normal distribution with mean and standard deviation equal to the observed X7 count-rate and its error, respectively. With this approach, for a given energy band, we created a simulated light curve with the same length as the original one, we performed a Lomb-Scargle analysis on it and we determined the main peaks of the periodogram. We simulated 10 000 light curves in the three energy bands and drew the distribution of the Lomb-Scargle peaks. For each peak, we estimated the mean and the standard deviation of its distribution by fitting it with a Gaussian model. The standard deviation was finally associated with the uncertainty on the periodicities determined from the observed light curves. The uncertainties (1σ) on the periodicity for the 0.3–10 keV, 1–10 keV and 0.3–1 keV are 2.4 d, 2.4 d and 3.4 d, respectively, which makes all the periodicities widely compatible with 240 d within less than 2σ. We used 240 d as reference time for the next analyses.

Should the shorter periodicity at ∼240 d be real in these data, we note that it is unlikely to be sinusoidal in shape. The agreement of the best-fit sinusoidal curve and the X7 light curve is not very suitable for the high flux states of the source and it tends to deviate at the epochs where the Lomb-Scargle fit was not performed. To confirm the non-sinusoidal behaviour, we run also an epoch folding search with the FTOOLS task EFSEARCH in the time range before January 2023 of the 1–10 keV light curve. We searched around a period of 240 d, binning the curve with 20 bins and with a time resolution of 10 000 s. Although the light curve has numerous time gaps which can introduce spurious signals, we found the strongest periodicity at 240.15 d. Fitting the χ² profile with a constant plus a gaussian, we estimated that the error on the period is 0.1 d at 1σ. This value is therefore compatible with the LS one. We folded all the Swift/XRT observations in the 0.3–10 keV band with the a 240 d periodicity (taken, as average value; see Fig. 2-bottom), indicating that the folded profile presents a peaked signal followed by a flux decay. In Fig. 1, we show the complete X7 long-term light curve where, in the inset, we overlaid the profile of the ∼240 d periodicity.

In order to test the statistical significance of the signal found with the LS, we firstly performed an approach based on scrambling the data of the light curve. Specifically, we assigned each observation a count rate taken from another, random observation; on the new generated curve, we performed a LS analysis and we took the maximum peak power in the period range of our signal. We created 10 000 new light curves for the 1–10 keV energy range (taken as reference) and we found that the number of signals with a peak power equal or higher than 0.229 (the power found in the real data for the periodicity) was < 10, implying that the ∼240 d observed signal is significant at > 3σ.

However, we note that the scrambling of the date and the false alarm probability of the LS do not take into account any red-noise, that can be important with long-term periodicities (Scargle 1981) and can lower the significance of the signals (Vaughan et al. 2003). We further created a set of 10 000 simulations following the approach presented in Sect. 3 of Walton et al. (2016). We simulated a light curve assuming that the power spectrum has a power-law shape with an index of one or two, with the index values being a guess on the real slope. The light curves were simulated with a time bin of 2 ks, i.e. the average exposure of the Swift/XRT monitoring. The length of the simulated light curve was increased to 20 times the observed light curve. Then, we picked only the time bins consistent with the times of the real Swift/XRT observations. For each filtered light curve, we performed a LS analysis. We then determined how many peaks have power equal to or larger than the one found in the observed data. We calculated it for the 1–10 keV energy band, as reference. We found that, between 200 d and 300 d, a total of 87 and 1042 peaks over 10 000 were found above a power of 0.229 for a power-law index equal to one and two, respectively. Hence, in the presence of red noise, our results would imply that the actual significance of the 240 d periodicity could be lower than the 3σ reported above.

These tests prove that the real nature of the periodicity cannot be verified yet and more observations are necessary.

3.2. Short-term temporal analysis

The XMM1 and XMM2 observations of X7 show that, in both epochs, the source was in a similar regime, with strong and repeated flares (Fig. 3). The flares enhance the observed emission of a factor of two with respect to the persistent level and, more interestingly, they all appear flat-topped at a very similar flux. We note that the duration of the flares is not regular and varies from a few to tens of kiloseconds.

Fig. 3. Background-subtracted light curve of X7 in XMM0, XMM1, and XMM2 (from left to right), binned at ΔT = 300 s (see the top-left legend), for the 0.3–10 keV (top panel), 0.3–0.8 keV (centre-top), and 0.8–10 keV (centre-bottom) energy bands. We also show the hardness ratio between the hard and the soft band (bottom panel). The dashed vertical lines divide the data between XMM0, XMM1, and XMM2, where the real time gaps between them were removed for displaying purposes only.

Time lags. We performed a search on time lags for XMM2, over the same timescales and energies as those reported for XMM1 in Pintore et al. (2021). However, no significant lags were found, implying that the hard lag detected in XMM1 is not observed again.

Search for quasi periodic oscillations. We searched for quasi-periodic signals (QPO) in XMM2⁴ data by creating power density spectra (PDS) from segments of 3000 s each, resulting in a frequency resolution of approximately 0.0003 Hz. We selected the length of these intervals to ensure sufficient photon accumulation for an adequate signal-to-noise ratio, while also preserving any potential transient signals. We produced a total of 27 PDS.

In some PDS, a narrow feature at about 1 mHz was observed, although it was not statistically significant. Similar features are sometimes visible in photon-starved EPIC-PN datasets, which may indicate that they are either statistical fluctuations or instrumental artefacts. We also combined the 27 PDS into a single averaged PDS to reduce uncertainties. This choice is supported by the fact that each PDS is very poor in statistics and there are no significant differences between the PDSs extracted in segments containing either the flares or the persistent emission. The PDSs are all dominated by noise, implying no significant variability over timescales ≤3 ks. The averaged spectrum is dominated by Poissonian noise as well, although an unresolved peaked feature is visible at roughly 0.5 mHz. This frequency is very close to our frequency resolution, so we conclude that it is unlikely to be a real signal.

Pulsation search. We then conducted a search for coherent signals from X7 in XMM2⁵.

We performed accelerated searches for periodic signals: to account for potential variations in spin period due to intrinsic spin-down or spin-up, the proper motion of the system, or the orbital motion of the compact object, we initially corrected the arrival times of the events using a grid of approximately 930 Ṗ/P values within the range of (±) 10⁻¹¹ s s⁻² to 10⁻⁵ s s⁻². We note that this does not take into account any orbital correction. The search did not yield any statistically significant signals and resulted in a 3σ upper limit of 5% on the fractional amplitude of a signal within the period range of 0.146 to 100 seconds.

Subsequently, we conducted an extensive search by applying binary-system orbit corrections to the photons ToA, accounting for the motion of the compact object around its companion. This involved the analysis of ∼30 000 different orbital configurations with orbital periods ranging from 4 hours to 4 days and orbital separations ranging from 2 to 120 light-seconds. A marginally significant signal was detected over several orbital corrections corresponding to a spin period P ≃ 2.73 s (Fig. 4). Its statistical significance could be reduced further due to the large number of trials. However, the large number of orbital parameters that gave rise to this signal detection encouraged us to search for such a signal in the XMM1 data. In fact, we did detect a quite close signal at P ≃ 2.63 s in this observation. Even though the detected signal in XMM1 was also marginally significant, the result that a very similar signal was detected in two separate observations improves its statistical significance. Taking into account that the signal has a very low pulsed fraction⁶ (PF) of ∼5%, it means that a low power in the PSD is expected. The orbital parameters are quite unconstrained, providing an orbital period in the range 2–4 d and A_x sin(i) of ∼20–120 lt-s. The detected signals would imply a quite strong, three-year spin-down rate of Ṗ/P ∼ −1.06 × 10⁻⁹ s s⁻¹. Still, the maximum significance of the signals is close to ∼3.5σ, which does not allow us to firmly claim it is a real signal. Therefore, we present it as a possible/candidate spin period. Further observations are needed to confirm its significance.

Fig. 4. Estimates of the X7 orbital period (Porb) and the projected semi-axis of the orbit for a spin period of ∼2.73s in XMM2. The confidence contour, more similar to an atoll rather than a ‘banana’ shape, is indicative of the low significance of the spin period detection. However, a similar plot can be obtained for XMM1 as well, boosting the real significance of the signal. The colour bar indicates the Leahy power of the signal.

Pulsation search with evolutionary algorithms. We also conducted another search for coherent signal in XMM1 and XMM2 using stingray (version 2.0.0) and particle swarm optimisation (PSO), an evolutionary algorithm (Kennedy & Eberhart 1995), similar to the search presented in Sacchi et al. (2024). PSO is a meta-heuristic algorithm designed to find the maximum value of a dependent variable in a multivariate function. In our case, the independent variables are the two orbital parameters (A_x sin(i) and orbital period), Ṗ/P, and the phase, while the dependent variable is the power of the signal. These independent variables define a ℝ⁴ space (the search space). More details about the method are described in Appendix A.

The parameter space for A_x sin(i) and the orbital period P_orb present an evident “banana” shape for two signals at ∼2.63s and ∼2.73s for XMM1 and XMM2, respectively, hence confirming the results obtained with more standard algorithms. These signals are associated with a maximum power of 47.99 for XMM1 and 45.27 for XMM2. We observed that high power values show a wide range of $\dot{P}/P$ ([−1 × 10⁻¹², 2 × 10⁻¹¹] s s⁻¹ and [−3 × 10⁻¹², 1 × 10⁻¹⁰] s s⁻¹ for XMM1 and XMM2, respectively), which implies that Ṗ cannot be constrained. On the other hand, we found that the uncertainties (at 1σ) for the coherent signals are P_spin = 2.628s ± 0.001s for XMM1 and P_spin = 2.726s ± 0.012s for XMM2.

It is important to note that the shapes of the two “bananas” obtained by these signals are totally overlapping (Fig. 5). Such a result suggests that the two signals may have a common origin since, if they were fake, it would be extremely unlikely to observe the same parameter space for the orbital parameters. To investigate if the coherent signal can be an artefact, we selected the largest number of background events (∼30 000) from XMM1 taken from the same CCD of the source and we performed again the analysis on these events. We did not find any signal at the same frequencies with powers comparable with (or higher than) those found for X7. The results support and strengthen the presence of pulsation in the source.

Fig. 5. Parameter space of the orbital parameters for the signal at 2.63s and 2.73s for XMM1 (left) and XMM2 (right), respectively. For displaying purposes only, in each plot the contours indicate the solution of the other observation. The colour bar indicates the Leahy power of the signal.

Fig. 6. Top: Hardness-intensity diagram for PN+MOS light curves of all XMM-Newton observations (after smoothing the data, see text). The data of the XMM0, XMM1 and XMM2 observations are reported in blue, orange and green, respectively. Hardness is defined as the ratio between the bands 2.0–10 keV (H) and 0.3–0.8 keV (S), while the rate is the sum of both bands. Smoothing with a moving window has been applied to the data (see text), and therefore we do not plot the errors on each point. The coloured and numbered segments on the top indicate the hardness ranges used for the spectra extraction shown in Fig. 8. Bottom: Hardness-intensity diagram for the available Swift/XRT data. Hardness is the ratio between the bands 1.0–10 (H) keV and 0.3–1.0 (S) keV.

Fig. 7. Combined EPIC-pn + MOS(1,2) spectra of the six spectral states extracted from the HID (see text) unfolded through a power law of index zero. The colour convention is the same as indicated in Fig. 6-top.

Fig. 8. Combined EPIC-pn + MOS(1,2) spectra of the six spectral states extracted from the HID (see text and Fig. 7). The colour convention is the same as indicated in Fig. 6-top. The best fit is given by the TBABS(DISKBB + DISKPBB + CUTOFFPL) model (see Table 2), where the DISKBB, DISKPBB and CUTOFFPL models are indicated with dashed, dot-dashed and dotted lines, respectively. In the bottom panel of each plot, it is shown the ratio of data and model.

3.3. Hardness-intensity analysis

We built a hardness-intensity diagram (HID) with all the XMM-Newton observations by using the EPIC data (pn + MOS1, 2) in the energy bands 0.3–0.8 keV and 2.0–10 keV (Fig. 6-top). The two bands are slightly different than those used in the previous section. With this choice, we could evaluate the HID using energies where we are confident the low and high energy spectral components of X7 dominate (i.e. below 1 keV and above 2 keV; Pintore et al. 2021) and therefore maximise the possibility of finding evidence of well-distinguished spectral states. The background-subtracted light curves, binned with ΔT = 50s, were smoothed with a moving window method in order to reduce the noise fluctuations of the data, following the same approach proposed in D’Aì et al. (2021). The HID shows the existence of at least three spectral regimes (Fig. 6-top), two sharing a similar hardness but different fluxes and one which is both harder and brighter than the others. We labelled the lowest luminosity state (which entirely represents the observation XMM0) as ‘low state’, and the ‘persistent state’ and ‘flaring state’ are the other two. The two latter states are populated in a very similar way in both XMM1 and XMM2. We also mention the existence of a branch connecting the ‘persistent state’ and the ‘flaring state’, which indicates a smooth transition between these two states (albeit a rapid one, given the low density of these data points).

We checked if this behaviour is also observable in the Swift/XRT data. Because of the low count-rate in Swift/XRT, we were forced to use the energy bands 0.3–1.0 keV and 1.0–10 keV to have enough counting statistics in both bands. Despite this precaution, the uncertainties due to the poor statistics are still large and do not allow us to robustly constrain the existence of the three regimes (Fig. 6-bottom); however, we note that a similarity between the Swift/XRT and XMM-Newton HIDs can be found, although we warn the reader that the baseline of all the Swift/XRT observations cover a significantly longer time range than the single XMM-Newton observations. Furthermore, we highlight that it might be possible that flares similar to those observed in XMM-Newton were caught by Swift/XRT only partially due to the short exposures (in general, around 2 ks), affecting the final shape of the HID. Finally, the different energy bands used here and the sensitivity of each instrument can also be responsible for a different HID in Swift/XRT.

3.3.1. Spectral analysis

We focus our spectral analysis on the high-quality XMM-Newton observations. The HID tracks very well the rapid spectral evolution of the source over time and allowed us to investigate in detail the flare peaks, the persistent emission, and the transition phase between these two regimes.

Because the X7 spectral properties in XMM0 were already reported in Pintore et al. (2021) and are not particularly informative, we focussed only on XMM1 and XMM2. We selected six time intervals associated with hardness < 0.4 (the persistent emission), [0.4–0.5], [0.5–0.8], [0.8–1.15], [1.15–1.5] (the transition) and > 1.5 (the flares), and we extracted the spectra for both XMM1 and XMM2. This reflects to a total of 12 EPIC-pn and 24 EPIC-MOS (1 and 2) spectra. After checking that similar spectral states have consistent pn and MOS spectra, we combined them with the SAS task EPICSPECCOMBINE, creating a final number of 6 EPIC (pn + MOS1, 2) spectra (Table 1 for details). We performed the fits with XSPEC v.12.13.0c (Arnaud 1996) in the energy range 0.3–10 keV. In Fig. 7, we plot the six spectra unfolded with a power-law of index zero in order to clearly show the source spectral evolution.

We described the data with a model already reported as suitable for X7 (Pintore et al. 2021), namely, the combination of two thermal components and a cut-off power law (CUTOFFPL in XSPEC), absorbed by a neutral column density (TBABS in XSPEC). The two components are a DISKBB and DISKPBB⁷ model (e.g. Mineshige et al. 1994; Kubota & Done 2004). The cut-off power law was found in a number of ULXs, and it might hint to an accretion column above a magnetised NS (although other scenarios are possible; e.g. Mills et al. 2024). We note that the additional cut-off power law was also found to be necessary in the spectra of NGC 4559 X7, when a combined XMM-Newton and NuSTAR spectral analysis was performed (Pintore et al. 2021).

3.3.2. Model 1: DISKBB+DISKPBB+CUTOFFPL

Since we found no significant variability between the individual spectra, we tied the column density (modelled with TBABS) and the p parameter of the DISKPBB amongst all the spectra. The CUTOFFPL parameters were instead frozen to Γ = 0.59 and E_cut = 7 keV (i.e. the average of the cut-off power-law parameters of the PULXs; Walton et al. 2018). We left free to vary all the remaining spectral parameters of the two thermal components and the normalisations of the CUTOFFPL model. The best-fit gives χ_ν² ∼ 1.06 for 3776 d.o.f. (null hypothesis probability (NHP) = 4 × 10⁻³) and we present the results in Table 2). Although it is formally not statistically acceptable, we note that the continuum is generally well modelled and no evident residuals pop up, with the exception of a few narrow residuals around 1–2 keV (see Fig. 8).

The disc temperature and its bolometric luminosities are often compared in order to put constraints on the nature and geometry of the accretion disc itself. Indeed, for a standard, thin accretion disc, the luminosity and the inner disc temperature are expected to correlate as L_x ∝ T⁴ (Shakura & Sunyaev 1973), for a constant emitting area (i.e. constant inner disc radius). Instead, in the case of an advection dominated disc, the correlation changes as L_x ∝ T² (Watarai et al. 2001). Therefore, looking for such correlations can give rough hints on the structure of the disc and its accretion regime. A number of works investigated the existence of these correlations in many ULXs (e.g. Feng & Kaaret 2006; Kajava & Poutanen 2009; Pintore & Zampieri 2012; Walton et al. 2013, 2020; Robba et al. 2021; Barra et al. 2024) which resulted in a variety of relations, although we note that the outcomes are often dependent on the continuum spectral modelling. However, the search for correlations between luminosity and disc temperature is nevertheless a useful tool to characterise the accretion flow. Hereafter, we look for correlations for the thermal components used to model the spectral continuum of NGC 4559 X7.

The soft thermal component slightly changes in flux but does not show significant variations of its temperature, apart for the highest flux state where such a component is very weak and poorly constrained (kT_cold ∼ 0.35 keV). We found a tentative $L_{\text{cold}}\sim {\text{kT}}_{\text{in}}^{-2}$ anti-correlation, but the uncertainties are large in both normalisation and index; therefore we caution that the anti-correlation may not be correct. The harder thermal component has a temperature higher than ∼1 keV and it shows instead a correlation that is compatible with $L_{\text{hot}}\sim {\text{kT}}_{\text{in}}^{1.37\pm 0.05}$ (Fig. 9-left). The apparent radii span the ranges ∼1000 km − 4500 km and ∼70 km − 250 km for the soft and hard thermal component, respectively (Fig. 9-right), assuming a low inclination of the system (10°). We found that the apparent radii are not constant in time but that can vary with the source flux. The L_x ∝ kT⁴ correlation is given assuming the disc at a fixed inner radius (e.g. Kubota & Done 2004; Done & Kubota 2006). Hence the correlation found between luminosity and temperature can be strongly affected by the radius change. We highlight that the XMM0 data, fitted with the same model, are located in the L_x–kT region which is unrelated to the others. It indicates that the flux changes are key to the source evolution.

Fig. 9. Thermal component unabsorbed bolometric luminosities (0.001–30 keV) versus the soft (blue points) and hard (orange points) disc temperature (left), and thermal component temperatures (soft disc in blue, hard disc in orange) versus their inner disc radii (right) for Model 1. We assumed an inclination of the system of 10 degrees. The circle, triangle, diamond, cross, star, and square points indicate the spectra from one to six indicated in Table 1, respectively.

Notably, the cut-off power law still dominates the very high energies and its emission correlates with the global flux. We estimated that the flux ratio between the hard components (hard thermal disc plus cut-off power law) and the soft thermal components increases from 1 to ∼6, moving through the lowest to the highest flux spectral state, respectively. In particular, we found that the soft component manifests a maximum variation of a factor of ∼1.5 (taking into account the uncertainties), compared with the factor of ∼3.5 of the hard ones, with the latter mainly being responsible for the spectral changes (as seen also from the light curves in Fig. 3). We further note that the DISKPBB luminosity and cut-off power-law luminosity are strongly correlated, suggesting that their spectral changes are triggered by the same physical process.

3.3.3. Model 2: DISKBB+DISKPBB+CUTOFFPL with a constant radius

The previous model gave a complex evolution of the spectral shape of X7 and indicated that the emitting radii of the thermal components change during time, although by no more than a factor of 2–3 taking into account the quite large uncertainties. Therefore, we finally tested a model where the thermal component normalisations of both the DISKBB and DISKPBB are tied between the spectra (i.e. assuming fixed emitting radii). The best fit (χ_ν² ∼ 1.1 for 3789 dof) is shown in Table 3. It converged towards radii of ∼3500 km and ∼50 km for the soft and hard component, respectively. Remarkably, the soft component shows a robust correlation between the inner temperature and its bolometric luminosity ($L_{\text{cold}}\sim {\text{kT}}_{\text{in}}^{2.9\pm 0.3}$; Fig. 10), which is neither compatible with a standard accretion disc nor with a fully advection dominated disc. We note that the temperature of the cold component is at the highest during the minima of the global source flux and the component is negligible during the peaks of the flares, i.e. the component anti-correlates with the X7 bolometric flux. On the contrary, the behaviour of the hard component is less constrained ($L_{\text{hot}}\sim {\text{kT}}_{\text{in}}^{4.6\pm 0.9}$; Fig. 10), and compatible with a standard disc although the p parameter of the DISKPBB indicates an advection dominated disc (p ∼ 0.5 − 0.6). These contradictory results might suggest that the assumption of a fixed inner radius might be suitable for the cold component but not for the hot one. We also highlight that the cut-off power-law component is negligible at the lowest X7 fluxes (≪1% of the total flux), while it is quite strong during the flares (∼25%).

Fig. 10. Thermal component unabsorbed bolometric luminosities (0.001–30 keV) versus the soft (blue points) and hard (orange points) disc temperature for Model 2. The circle, triangle, diamond, cross, star and square points indicate the spectra from one to six indicated in Table 1, respectively.

For completeness of our spectral analysis and in addition to the models described above, we also tested an absorbed DISKBB+DISKPBB, DISKPBB+DISKPBB+CUTOFFPL and a DISKPBB+CUTOFFPL, but they gave poorer or less constrained results. For this reason, they are not taken into account in this work.

3.4. High-resolution spectroscopy

We found that some residuals appear around 1 keV in the EPIC spectra, especially for the highest flux state (purple and cyan spectra in Fig. 8) and more evident in spectra extracted from the entire observation (see Pintore et al. 2021 for the archival data). Such features are common in ULX CCD spectra taken with Chandra and XMM-Newton satellites (e.g. Middleton et al. 2014). The combined RGS spectra along with the corresponding EPIC-pn and MOS 1, 2 spectra are shown in Fig. 11. The broad feature seen in the EPIC data around 1 keV is resolved into a few emission lines in the RGS spectra, which provide a high resolving power R_RGS = E/ΔE ∼ 200 (400) at 1 keV for the 1^st (2^nd) spectral order. Further, weaker, emission features seem to be present. We focussed on the RGS data obtained by stacking the two longest observations (XMM1 and XMM2 as mentioned in Sect. 2). Besides, XMM0 has a different and fainter spectrum which makes it more prone to issues with the background subtraction.

Fig. 11. Combined RGS first and second order (grey and blue), EPIC-pn and MOS 1,2 (black, red, and green) spectra (top panel). The black and red lines refer to the continuum and continuum plus CIE models, respectively. The corresponding residuals are shown in the middle and bottom panels. We note that the Ne K – Fe L emission complex is near 1 keV.

Following previous works (see e.g. Pinto et al. 2021), we kept the EPIC data only between 2–10 keV, i.e. outside the RGS energy, owing to the much poorer resolution (R_EPIC = E/ΔE ≲ 20 below 2 keV, which is not sufficient to resolve narrow lines). This avoids introducing a further degeneracy between different models that would provide comparable spectral fits of the EPIC data (its count rate higher than the RGS would basically lead the fit). The energy range between 0.5–2 keV (0.8–1.8 keV) was chosen for the RGS first (second) order because outside this range the spectrum is background dominated.

The RGS+EPIC spectral fits were performed with the SPEX fitting package v3.07 (Kaastra et al. 1996), which is equipped with up-to-date atomic databases and flexible plasma models. The spectra were grouped according to the optimal binning discussed in Kaastra & Bleeker (2016) and the Cash statistics (C-stat; Cash 1979) was used. This has the advantage of smoothing the background spectra in the energy range with low statistics and removing narrow spurious features. We note, however, that further tests that were performed by rebinning to at least 25 counts per bin provided consistent results.

Following the modelling of the HID-selected EPIC spectra, we adopted a broadband continuum model consisting of two thermal components (disc-blackbody, DBB⁸ in SPEX.) and a power law with a high-energy cut-off (i.e. a POW component multiplied for an ETAU component with Γ_pow = 0.59 and E_cut = 7 keV, respectively). The continuum model was self-consistently fit to the five averaged RGS+EPIC spectra except for a multiplicative constant that was kept free to vary with respect to the EPIC-pn (which was fixed to 1) to account for the typical 5% cross-calibration uncertainties. The RGS+EPIC continuum model provides an overall C-stat of 811 (and a comparable value of χ²) for 659 d.o.f. The continuum model is represented by the solid black line in Fig. 11 with the corresponding residuals in the middle panel. The precise form of the continuum model has little impact on the modelling of narrow emission feature (at least in the case of collisionally ionised plasma; see below).

3.4.1. Gaussian line scan

Some narrow features and residuals can be seen in Fig. 11 although they are clearly weaker than in other ULXs exhibiting softer X-ray spectra (e.g. Pinto et al. 2017); this is not surprising given that, in hard ultraluminous X-ray sources (Sutton et al. 2013), there is an overwhelming continuum from the innermost regions (see e.g. Pinto et al. 2017; Kosec et al. 2021). For a more efficient visualisation and a clearer detection of narrow lines, we performed a standard scan of the spectra with a moving Gaussian line following the approach used in Pinto et al. (2016). We tested five logarithmic grids with energies between 0.5 and 10 keV and different sizes according to the adopted line widths (300, 500, 1000, 1200 and 2000 points for a line FWHM of 10 000, 5000, 2500, 1000 and 100 km s⁻¹, respectively). At each energy, the ΔC improvement to the best-fit continuum model is recorded. Its square root provides the single-trial significance; to distinguish between emission and absorption lines, we multiplied the $\sqrt{\mathrm{\Delta }C}$ by the sign of the Gaussian normalisation.

The results of the line scan obtained for the combined RGS+EPIC spectra are shown in Fig. 12. We remind the reader that in the 0.5 − 2 keV energy band, only RGS data is considered. The line scan picked out the strong emission feature near 1.0 keV, which is broadly consistent with the Ne X rest-frame laboratory energy and some highly ionised iron species (Fe XX-XXIV). Further, weaker, emission features are close to the O VIII (0.65 keV) and Si XIV (2.01 keV) lines. The absence of lower-ionisation ionic species such as Ne IX and O VII, frequently found in the RGS spectra of softer ULXs, would suggest that the plasma is either highly ionised or mainly in collisional-ionisation equilibrium (CIE). Some absorption features, albeit weaker, appear between 0.6 − 0.8 keV, which mimic those previously found in many ULXs (see Sect. 3.4.2 and 4 for more detail). The RGS line search for the broad 10 000 km s⁻¹ exhibit an overall shape which mimics that of the EPIC residuals (see Fig. 8).

Fig. 12. Gaussian line scan for the combined RGS (0.4–2 keV) and EPIC (2–10 keV) spectra and different line widths. Rest-frame centroids of the typical dominant H-like and He-like transitions are labelled.

3.4.2. Physical models scan

In order to boost the significance of the detection of the plasma component(s) by simultaneously accounting for multiple lines, we performed a search throughout a multidimensional parameter space of plasma in photo- (PIE) or collisional (CIE) ionisation equilibrium. This technique prevents the fits from getting stuck in local minima, although it is computationally expensive (see e.g. Kosec et al. 2018b; Pinto et al. 2020). Given the weakness of the absorption features we only tested PIE absorption models as the results would not be distinguishable from CIE absorption models. Besides, absorption lines from CIE plasmas are not common in XRBs (see e.g. Pinto et al. 2021; Neilsen & Degenaar 2023). We adopted solar abundances to save computation time; this might have had an impact on the overall spectral improvement as the abundance pattern could differ from the solar one (see e.g. Barra et al. 2024 for the relevant case of ULX Holmberg II X-1).

CIE line emission – We first focussed on the emission features by running a multidimensional automated scan with an emission model that assumes collisional ionisation equilibrium (CIE model in SPEX). The CIE was multiplied by a redshift (REDS model) with FLAG=1 in order to test various line-of-sight (LOS) velocities (v_LOS). We adopted a logarithmic grid of temperatures between 0.1 and 5 keV with 50 points. The range of v_LOS was chosen between −0.3c (i.e. blueshift or outflowing plasma) and +0.3c (redshift or inflowing plasma). The steps in velocities were selected according to the adopted velocity dispersion or line broadening (v_step = 500, 750, 1000, 1500 and 3000 km s⁻¹ for v_RMS = 100, 1000, 2500, 5000 and 10 000 km s⁻¹, respectively). The CIE normalisation or emission measure was the only additional free parameter. Similarly to the gaussian scan, the ΔC improvement to the best-fit continuum model is recorded for each grid point and plotted in a kT − v_LOS graphic. The best results (ΔC_max = 21) are obtained for v_RMS ranging from 100 to 1000 km s⁻¹ with the case for 100 km s⁻¹ shown in Fig. 13 (top panel). The peak or best fit is achieved at a temperature of 1.81 keV and a v_LOS = 0.0253 c or 7600 km s⁻¹ (i.e. weakly redshifted) although a secondary, less preferred, solution is found around −0.015 c. In Fig. 11, we also show the best-fit model, including the CIE component (solid red line in the top panel and corresponding residuals in the bottom panel). This component mainly accounts for the Ne K and Ne L lines around 1 keV. The CIE has a X-ray luminosity L_{[0.3 − 10 keV]} ∼ 1 × 10³⁹ erg/s (for the adopted distance of 7.5 Mpc). The detection significance has to account for the look-elsewhere effect: in previous works, for RGS spectra with comparable statistics and with an identical parameter space explored, it was shown that a ΔC_max = 21 corresponds to a conservative confidence level of about 3σ (see e.g. Pinto et al. 2017, 2020; Kosec et al. 2018a,b). This has been further confirmed with an alternative method which uses cross-correlation instead of direct spectral fitting (see e.g. Kosec et al. 2021; Xu et al. 2022).

Fig. 13. Top panel: Multi-dimensional grids of emission models of plasma in collisionally ionisation equilibrium for the combined RGS (0.4–2 keV) and EPIC (2–10 keV) spectra. Middle and bottom panels: Grids of models of photoionised plasma in emission and absorption, respectively. Both scans adopted a line broadening of 100 km s−1 and solar abundances. The best-fit CIE model is also shown in Fig. 11.

PIE line emission – We also attempted to describe the emission lines with models of plasmas in photoionisation equilibrium. SPEX is provided with the unique PION model, which instantaneously calculates the photoionisation balance of the gas by using the fitted continuum model, produces the emission lines and recombination continua, and fits them to the spectra. Here, the balance is parameterised by the ionisation parameter, ξ = L_ion/(n_H R²), where L_ion is the ionising luminosity (typically measured between 13.6 eV and 13.6 keV or 1–1000 Rydberg), n_H is the hydrogen number density and R is the distance from the ionising source. In analogy with the CIE, we created grids of PION emission models adopting a logarithmic grid of ionisation parameters (log ξ between 0 and 6 with steps of 0.1). For the v_LOS and v_RMS we used the same ranges tested over the CIE grids. The only free parameter for the PION component is the column density, N_H. To produce only emission lines, we chose a unitary solid angle Ω = 1 and a null covering fraction f_cov = 0. The PION search resulted in a much lower improvement with respect to the CIE grids (ΔC_max = 11 instead of 21) but also shows a peak at v_LOS = 0.0253 c with log ξ = 2.9 (see Fig. 13, middle panel).

PIE line absorption – Finally, we also performed a scanning of grids through absorption models by selecting a solid angle Ω = 0 and a covering fraction f_cov = 1 for the PION component. Since we do not expect redshifted absorption for the winds, we only selected a range of blueshifts for the LOS velocity, again following previous work. PION in absorption works in the same way as in emission; the main difference is that the emission lines component is additive, while the absorption line component is multiplied for the continuum components. As for the emission, the only additional free parameter for each grid fit is the PIONN_H. The absorption model grids provided worse results than the CIE grids (see Fig. 13, bottom panel), which is expected given the weakness of the individual absorption features (see Fig. 12). The best fit is achieved for v_LOS = −0.06 c and log ξ = 1, which would suggest a weak outflow, but the detection is not very significant given the maximum fit improvement of ΔC_max of just 14.

4. Discussion

NGC 4559 X7 is a ULX with structured temporal and spectral variability in the form of random flares, which appear to be all flat-topped, and a marked spectral variability characterised by “islands” in the HID. These properties, not very common in the ULX population, make NGC 4559 X7 an intriguing case of study.

4.1. Accretion rate variability or accretion geometry changes

During the almost 20 years of X-ray observations, the source showed a significant and prolonged variability on both short (hours to days) and long-term (months to years) timescales. The long-term variability is characterised by a flux evolution of at least a factor of 5, with rapid increases and subsequent decays on timescales of a few weeks. We estimated that the peak X-ray luminosity is ∼4 × 10⁴⁰ erg s⁻¹, while the average one is ∼5 × 10³⁹ erg s⁻¹, for a distance of 7.5 Mpc. The nature of such a variability might be explained with changes in the accretion rate or in the accretion geometry, or a combination of the two.

The source emission shows at least three spectral states, two characterised by similar shapes but different fluxes (that we labelled as ‘persistent’ states) and one (the ‘flaring’ state) very hard and brighter than the others. There is now an increasing group of ULXs which present these ‘island’ spectral states. Of this group, we mention the source M81 X-6, which was observed mainly in two spectral states with a characteristic harder-when-brighter behaviour (Weng & Feng 2018; Gúrpide et al. 2021) and it was suggested to host a weakly magnetised NS (Amato et al. 2023). Recently, a long-term spectral evolution was investigated for Holmberg II X-1 (Gúrpide et al. 2021; Barra et al. 2024) and it was reported that the source can switch between hard and soft states directly related to the source luminosity. A possible explanation could be linked to a changing view of the inner regions of the accretion flow, shielded from time to time by an optically thick outflow ejected by the accretion disc. This is regulated by the closing of the funnel created by the outflow for increments in the accretion rate. For Holmberg II X-1, a super-soft state was also observed, while for NGC 4559 X7 such a spectral regime was never seen. On the other hand, we cannot exclude that the sources contain the same compact object and slightly different accretion rates or inclination angle of the binary system.

In all the three states of NGC 4559 X7 the soft, thermal component does not change significantly in temperature but only in normalisation, of a factor of ∼2–2.5 (this work and Fig. 9 in Pintore et al. 2021). This suggests little accretion rate variability at large radii from the compact object. Instead, the hard component is very variable (more than a factor of 5) and it is responsible for the ‘harder-when-brighter’ behaviour in X7, indicating that either the accretion onto the compact object is variable or other processes affect its emission.

We can hypothesise that the global variability may be a combination of small fluctuations in the accretion rate and changes in the accretion geometry that outshines the inner regions emission, or as occultation from an outer, optically thick wind. Some evidence for winds in this source were found in the XMM-Newton data. Both EPIC and RGS (XMM1 and XMM2) show spectral features around 1 keV which in RGS are resolved into a few Ne K and Fe L emission lines. A few absorption features may also be present, although the significance of them is not particularly high. The physical modelling favours collisional ionisation or shocks for the emission lines with a low velocity (redshift) of a few thousand km s⁻¹, which basically means plasma close to rest, given the limited resolving power and uncertainties around 1 keV. This is in line with the expectations for an intermediate-hard spectrum ULX where the binary system is expected to be seen close to face-on (Middleton et al. 2015). In such a condition, winds are less prominent, especially in the case of an NS with high magnetic field as compact object, which would inhibit the launch of powerful, fast winds if the disc is truncated above the spherisation radius. However, in most cases absorption lines were detected in ULX spectra with a number of counts higher than those collected for X7. The line emission component has a X-ray luminosity L_{[LE, 0.3 − 10 keV]} ∼ 1 × 10³⁹ erg s⁻¹, which is a factor of 5–10 higher than in typical ULXs (Pinto & Walton 2023). Instead, the overall contribution to the source X-ray luminosity is a few percent, which is in line with other ULXs. This might suggest that either the adopted distance of 7.5 Mpc is an overestimation (quite unlikely because it is well constrained) or that the underlying emission process is particularly efficient in this source. The CIE best-fit solution indicates a redshift of 0.0253c or 7600 km s⁻¹ (very similar to some observations of the precessing jet in SS 433 – Marshall et al. 2013; Medvedev et al. 2019 – and the outflow in NGC 247 ULX-1 – Pinto et al. 2021). If the compact object was an NS, it might imply emission from near the magnetosphere or leakage towards the NS where gravitational redshift starts to kick in (rather than from a shocked disc wind). This might explain the unusual luminosity of the line-emitting plasma. Deeper observations are needed to place stronger constraints on the X7 features.

4.2. A tentative super-orbital variability

A tentative periodicity of ∼240 d was found in the Swift/XRT monitoring for the years ∼2020–2022. If it was orbital, the system would be quite large, making hard to explain the extreme mass transfer rate that sustains the accretion onto the compact object (unless it was an NS). A Be-binary system, with high inclination and eccentricity, could also explain the flux peaks (entrance in the Be-disc, i.e. higher accretion) and the lower flux state (the exit from the disc, i.e. lower accretion). A more likely scenario is that the periodicity is super-orbital as often seen in long-term monitored ULXs (e.g. Fürst et al. 2017; Walton et al. 2016; Brightman et al. 2019; Vasilopoulos et al. 2020; Salvaggio et al. 2022).

However, we note that the periodicity was not confirmed in the years 2023–2024 and that red-noise could also affect the significance of the signal. Therefore, we cannot robustly confirm yet the existence of a super-orbital variability in X7. On the other hand, the long-term light curve has large gaps where the behaviour of the source was unknown and we cannot rule out that X7 experienced abrupt flux decreases/increases during them. This may be important as sources with clear super-orbital periods (for example the PULXs NGC 5907 X-1 and NGC 7793 P13; e.g. Fürst et al. 2021, 2023) did not recover immediately the periodicity right after a flux switch-off. Furthermore, it is not possible to exclude that the periodicity changes during time, as seen also in the Galactic source Her X-1 (e.g. Leahy & Igna 2010) or in the PULX M51 ULX-7 (Brightman et al. 2020, 2022).

4.3. The flaring activity

On short timescales, NGC 4559 X7 experienced a flaring-like activity, observed only at the very highest fluxes, i.e. in two XMM-Newton and one NuSTAR observations (Pintore et al. 2021). It cannot be proven yet if such an activity is directly connected to the high flux states, as long and continuous exposure observations are missing at low fluxes. The Swift/XRT monitoring did not catch any significant evidence of a flare, although we remark that the Swift/XRT exposures are short (on average, ∼2 ks) and the flares present a variable duration. The rising and decay times are quite fast (about 1 ks) and the peaks can last from a few hundreds of seconds up to tens of kiloseconds. This makes the flaring activity unpredictable and it is more probable that the Swift/XRT monitoring can observe only either the peak or the low flux state rather than the rising and decay epochs of the flares. Indeed, we estimate that the time spent during the rising/decay is on average ∼10% of the two high flux XMM-Newton observations (total exposure of 195 ks). However, we note that the HID of the Swift/XRT data shows a bunch of observations with (tentative, because of the large uncertainties) hard spectra at the lowest fluxes. Hence, although the current evidence is quite weak, it cannot be ruled out that flares also happen at low fluxes.

A quite impressive behaviour of the flares is that they are all flat-topped at the peaks, suggesting that there is a common mechanism able to limit the emission above a certain luminosity. A clear spectral variability is observed during the flares, mainly driven by the hard emission (> 1 keV). The spectral analysis of X7 revealed that each spectral state can be described with the combination of two thermal components (disc-like emission) plus a non-thermal one (a cut-off power law). As mentioned before, the softer, thermal component is marginally varying during time between the two highest flux states (no more than a factor of ∼2 in flux), suggesting that there are no robust changes in the accretion rate (or at least not on the observational timescales) at large distance from the compact object. The soft component is quite sensitive to the continuum adopted to model the spectra, and it might show, alternatively, a correlation or an anti-correlation between temperature and luminosity. We suggest that the outer accretion flow could be a thick disc distorted by the outflow or the standard, outer disc itself.

The apparent size of such a component is quite similar to the spherisation radius of a 10 M_⊙ BH accreting at ∼10 times the Eddington rate. Indeed, the spherisation radius is defined as

$$\begin{array}{c}\hfill R_{\mathit{sph}}=27/4·\dot{m}·R_g,\end{array}$$(1)

and for R_g = 15 km and ṁ = 10, the radius is ∼10³ km. This value is even smaller if the compact object was an NS (R_sph ∼ 1400 km for ṁ = 100). We should keep in mind that a flux variability similar to that of X7 is observed during some regimes (i.e. the ρ class) of the BH GRS 1915+105 (e.g. Belloni et al. 2000). Therefore, we cannot exclude that X7 hosts a massive BH, in a scaled version of Galactic BH accretion regimes.

The harder components are the most intriguing spectral features of X7. They become more and more dominant while the flux increases (about a factor of 5 in flux variability). We note that the X7 spectra pivot around 0.8–0.9 keV (i.e. the luminosity of the soft component is anti-correlated with that of the hard component), indicating that the spectral evolution cannot be explained in terms of a dipping activity but rather to, again, either a change of the accretion geometry of the flow in the innermost regions or a switch-off and switch-on only of the accretion rate on the compact object. Precession effects from the wind could explain the spectral variability (e.g. Middleton et al. 2019), but it would be hard to reconcile it with the absence of a quasi-periodic variability of the flares. The thermal, harder component shows a correlation, irrespective of the adopted continuum, between temperature and luminosity ($L_{\text{hot}}\propto k{T}_{\text{in}}^{1-4.5}$). This may suggest that the inner flow is compatible with an accretion disc, whose proper nature (standard or advection dominated) cannot be robustly constrained, although the spectral fits support likely an advection dominated disc (p ∼ 0.6). The smallest apparent emitting radius of the hot thermal component was estimated around 70 km, during the brightest spectral state. If it was the innermost stable orbit of a BH, this would imply a M_BH ∼ 8 M_⊙ (assuming a non-rotating BH and an disc inclination of 10°). In order to give a better estimate, we take into account a colour-correction factor of f_col = 1.7 (e.g. Shimura & Takahara 1995) and a boundary condition factor ξ = 0.4 (e.g. Kubota et al. 1998) for the accretion disc: this would imply a larger apparent radius and the BH mass estimate will re-scale to ∼9 M_⊙. However, the assumption that links the apparent radius of the high energy thermal component to the innermost stable orbit could not be appropriate; therefore we still cannot rule out the presence of an NS. We remark that the spectral modelling we adopted for NGC 4559 X7 is consistent with that used for the PULXs (e.g. Walton et al. 2018).

4.4. A possible new candidate PULX

It is worth to note that we found some evidence for a pulsed signal with a period of P ∼ 2.6s − 2.7s. Although marginally significant in the single observations, it is noteworthy that similar signals were spotted in two different XMM-Newton observations (the ones at the highest flux). The signals present the same parameter space for the orbital parameters and, in addition, our tests exclude that they are an artefact. These findings strengthen the existence of the pulsation in X7. Applying orbital corrections, we estimated that the orbital period falls in the range 2–4 days with a quite unconstrained A_x sin(i) (∼20–120 lt-s) that, if confirmed, could support a super-orbital origin for the 240 d periodicity in X7. If the pulsed signal is real, the compact object would be unequivocally an NS accreting well above its Eddington limit (about a factor of 100). The orbital parameters are compatible with those observed in the PULXs (e.g. Israel et al. 2017a; Fürst et al. 2017; Rodríguez Castillo et al. 2020; Bachetti et al. 2022; Belfiore et al. 2024). On the other hand, we remark that our estimate of the significance may be affected by caveats; therefore, we maintain caution regarding claiming the discovery of a new PULX before further confirmation. However, from a spectral point of view, PULXs are generally very hard (e.g. Pintore et al. 2016; Gúrpide et al. 2021) similarly to what spectrally found for NGC 4559 X7 in XMM1 and XMM2, where the source might show evidence for pulsations. Hence, there are many indications that point towards an NS in the source.

We found that the spin period would be affected by a strong secular spin-down of Ṗ ∼ 10⁻⁹ s s⁻¹, which is very high and quite unexpected since the source was always observed in accretion. If the propeller⁹ was active in between our observing epochs, we should have detected a luminosity jump of at least one order of magnitude for the observed spin period of 2.6 s, which is not present in the long-term light curve though. Indeed, the maximum jump is of about a factor of five. However, it was proposed that the lack of advection in the accretion flow and a combination of magnetic dipole field and ṁ could imply small drops in luminosity during propeller (Middleton et al. 2023). Furthermore, it cannot be forgotten that numerous gaps are present in the long-term light curve, and even though there are no real clues about any switch-off of the source during them, we cannot rule out such a hypothesis.

A first-guess possibility that we might consider is that the largest part of the source flux does not come actually from the accreting matter falling into the compact object, but rather from the one beyond the magnetospheric radius. Indeed, if the source was continuously in a propeller phase (and the feature around 1 keV in the highest flux spectrum can be perhaps due to shocked gas produced by this effect), the vast majority of the emission would arise from the matter outside the magnetospheric radius. From time to time, there could be a leakage of material to the NS surface which is responsible of the hardening of the spectral emission, an increase of the flux during the flares and possibly the appearance of the pulsation. In such a scenario, the pulsation is expected to be intermittent but this cannot be proven with the current counting statistics. The accretion process might proceed through a different mechanism as well.

We estimate that the co-rotation radius for an NS spinning at 2.7s is R_cor = (GM_NSP²/(4π²))^1/3 ∼ 3300 km (where G is the gravitational constant and M_NS is the NS mass). This is roughly compatible with the spherisation radius that we estimated above, which would explain the weak winds due to an intrinsic low efficiency of the radiation pressure (Mushtukov et al. 2019). In such a condition, the hard components (the thermal disc and the cut-off power law) may be the emission from an accretion curtain inside the magnetosphere and an accretion column above the NS.

If the source was at the spin equilibrium, the NS dipolar magnetic field would be very high (∼10¹³ G), truncating the disc likely before the spherisation radius. However, should massive outflows be confirmed in the future, this would likely imply in turn that the source has a smaller magnetic field and its accretion proceeds in absence of propeller. The secular spin-down is therefore difficult to explain unless we take into account physical mechanisms such as the threading of the magnetic field at large radii, which exerts a net torque on the NS and is more effective than the spin-up torque due to the accretion. Recently, Bachetti et al. (2022) and Liu (2024) reported that the PULX M82 X-2 experienced periods of spin-down while accreting, suggesting that the source is at the spin-equilibrium and small changes in the accretion rate can trigger the propeller. Although there is no indication of a spin-equilibrium scenario for X7 (see above), the spin-down in NGC 4559 X7 might be similar to M82 X-2. Another pulsar with a strong secular spin-down while accreting is the Galactic GX 1+4 (see e.g. González-Galán et al. 2012, and reference therein), even though it is a low-mass X-ray binary probably unlike NGC 4559 X7 (Soria et al. 2004) and the other PULXs. It is noteworthy to remark that the spin-down of GX 1+4 is even stronger (∼10⁻⁷ s s⁻¹) than that of NGC 4559 X7. However, given the longer spin period of GX 1-4 than X7, the Ṗ/P for the two sources differs by no more than a factor of 3, which is quite impressive seen the profound differences in their properties. 4U 1538-522, a Galactic high-mass X-ray binary pulsar, shows also a clear spin-down without any flux variability (e.g. Hu et al. 2024, and references therein). LMC X-4 is also characterised by almost quasi-periodic variations of spin-up and spin-down episodes (e.g. Molkov et al. 2017).

Some models proposed that the geometry of the interaction between magnetosphere and accretion disc has an elliptical shape (Perna et al. 2006), suggesting that spin-up and spin-down can act simultaneously in different positions of the magnetosphere or that a retrograde accretion disc could be responsible for the pulsar spin-down during accretion (e.g. Makishima et al. 1988; Dotani et al. 1989; Nelson et al. 1997). However, the latter scenario is hard to reconcile with NGC 4559 X7, as how a retrograde disc can be sustained for many years would have to be explained. Therefore, no clear explanations could be currently proposed to describe the source behaviour.

4.5. An exotic scenario: A millisecond pulsar in NGC 4559 X7

If the reported spin period for NGC 4559 X7 is merely a statistical fluctuation, we propose a further, exotic scenario for this source. The flat-topped flares suggest that a common mechanism inhibits the accretion above a certain level, decaying rapidly towards lower fluxes with softer spectral shapes. If the rapid passage from the flares to a persistent emission was due to a propeller switch-on/switch-off, the luminosity jump ($\mathrm{\Delta }L\propto 170P^{2/3}{M}_{1.4M_{\text{sun}}}^{1/3}{R}_{{10}^{6}cm}^{-1}$; e.g. Campana et al. 2018) of a factor of ∼3 (or 5 if we consider the lowest source flux ever observed) would imply an NS spinning at a few millisecond periods (P = 2ms − 5ms). The dipolar Larmor luminosity (see Burderi & Di Salvo 2013, for a review) is

$$\begin{array}{c}\hfill L_{\mathrm{Lar}}=3.85\times {10}^{35}{B}_8^2{R}_6^6{P}_{-3}^{-4}\text{erg}{\text{s}}^{-1},\end{array}$$(2)

where B₈ is the dipolar magnetic field in units of 10⁸ G, R₆ is the NS radius in units of 10⁶ cm and P₋₃ is the spin period in units of 1 ms. This luminosity has a strong dependence on the spin period and the magnetic field. For a ms pulsar of a period in the range P = 2 ms − 10 ms, in order to have a luminosity significantly smaller than the one observed in NGC 4559 X7 (minimum luminosity of ∼5 × 10³⁹ erg s⁻¹) so that its contribution can be ignored, the dipolar magnetic field has to be < 10¹² G. However, if the source is seen at small angles, i.e. face-on, and the walls of the outflows can beam the inner region emission, the isotropic source luminosity would be smaller (for example, a beaming factor of 0.1 implies an order of magnitude less), and therefore the magnetic field reduces in consequence.

The co-rotation radius, defined as

$$\begin{array}{c}\hfill R_{\mathrm{cor}}={\left(\frac{G{M}_{\mathit{NS}}{P}^2}{4{\pi }^2}\right)}^{1/3},\end{array}$$(3)

for a ms pulsar of a 2–10 ms period, is in the range 20–80 km, which is compatible with the apparent emitting radius of the hot accretion disc found from the spectral analysis (taking into account also the colour-correction and boundary condition factors). In the case of spin-equilibrium, to explain the switch-on/switch-off of the propeller, the magnetic field would be of the order of 10¹² − 10¹³ G, to have a luminosity limit (e.g. Campana et al. 2018) compatible with the X7 luminosity peak (4 × 10⁴⁰ erg s⁻¹). Taking into account the beaming, the magnetic field will be lower than that.

Currently, this exotic scenario cannot be proven because there are no available X-ray observations with a ms time-resolution and a large number of collected counts to measure very small pulse fractions. We suggest that XMM-Newton observations in TIMING mode of this and many more bright ULXs without pulsation detections will be necessary to find evidence of such a short NS period.

5. Conclusions

We investigated the spectral and temporal properties of the ULX NGC 4559 X7. Variability up to a factor of approximately five on long-term timescales was observed thanks to a Swift/XRT monitoring. Also XMM-Newton caught a rapid short-term variability twice, whose behaviour resembles a flaring activity. Flares are mainly driven by changes in the hard emission (> 1 keV), while the softer emission is more stable in time (less than a factor of 2–2.5 in variability).

We analysed X7 spectra by adopting models based on the combination of two thermal components plus a cut-off power law, which is a spectral modelling often used for the PULXs. We found that such a continuum is suitable for the X7 spectra. The soft component is quite sensitive to the adopted continuum, and it is not clear if its temperature and luminosity follow a correlation or an anti-correlation. On the other hand, the hard thermal component correlates in the two quantities, irrespective of the adopted continuum. Hence, we suggest that the soft component is associated with the outer regions of the disc or a weak outflow above it, while the hard thermal component is likely an advection dominated accretion disc. The apparent emitting radii of the two components are ∼2000 km and ∼100 km, respectively.

A pulsed signal at 2.6s–2.7s was found in the XMM-Newton observations taken in 2019 and 2022. Although the signal is found in two observations separated by years and sharing the same parameter space for the orbital properties, on one hand it increases its significance, and on the other, the spin period variability implies a strong spin-down of ∼10⁻⁹ s s⁻¹. The latter is hard to explain, as X7 was always observed in an apparent accretion regime. Several hypotheses were taken into account, such as threading of the disc, a retrograde accretion disc, and a high NS dipolar magnetic field, but none appear to be conclusive. New XMM-Newton observations will be needed to definitively constrain the existence of such a pulsating signal.

Data availability

All of the data underlying this article are already publicly available from ESA’s XMM-Newton Science Archive (https://www.cosmos.esa.int/web/xmm-newton/xsa), NASA’s HEASARC archive (https://heasarc.gsfc.nasa.gov/), and the Chandra Data Archive (https://cxc.harvard.edu/cda/). The SPEX code is also publicly available at https://github.com/ciropinto1982.

Acknowledgments

This work has been partially supported by the ASI-INAF program I/004/11/4 and program I/004/11/6. FP, GR, CP, MDS, AD, AW acknowledge support from INAF grants OBIWAN (Observing high B-fIeld Whispers from Accreting Neutron stars, 1.05.23.05.12) and BLOSSOM. CP acknowledges support from PRIN MUR 2022 SEAWIND 2022Y2T94C funded by NextGenerationEU. MDS acknowledges support from the INAF grant ACE-BANANA (1.05.12.05.17). DJW acknowledges support from the Science and Technology Facilities Council (STFC; grant code ST/Y001060/1). NOPV acknowledges support from INAF and CINECA for granting 125 000 core-hours on the Leonardo supercomputer to carry out the project “Finding Pulsations with Evolutionary Algorithms”.

References

Appendix A: Particle swarm optimisation algorithm

The PSO algorithm searches for an optimal set of orbital parameters that yield the highest signal power in this space, without relying on a pre-determined grid. It relies on the movement of ‘particles’ (i.e. candidate solution in the search space), whose positions are updated iteratively. The position update relies on two formulas. The first defines how the velocity of a particle is updated:

$$\begin{array}{c}\hfill v(t+1)=\omega v(t)+c_p{r}_p(p-x(t))+c_g{r}_g(g-x(t)).\end{array}$$(A.1)

The second defines how its position is updated based on the velocity:

$$\begin{array}{c}\hfill x(t+1)=x(t)+v(t).\end{array}$$(A.2)

Here, ω is the inertia and controls the influence of the previous velocity v(t) on the new velocity v(t + 1), x is the particle position, p is the best position reached by the current particle (i.e. the position with the highest power), g is the best position reached by any particle, r_p and r_g are two random numbers, and c_p and c_g are acceleration coefficients and represent, respectively, the cognitive component (i.e. attraction towards p) and the social component (i.e. attraction towards g).

The algorithm consists of four steps:

• Initialisation: n particles are randomly distributed in the search space with a random initial velocity

• Evaluation: the power associated with each particle is computed

• While a stopping criterion is not met:

  • Update: the particles adjust their velocities and positions using Eqs. A.1 and A.2

  • Evaluation of the particles at the new positions

In this work, the search space was divided into six non-overlapping regions:

1. A_x sin(i)∈[20, 40] and P_orb ∈ [2, 3]

2. A_x sin(i)∈[40, 60] and P_orb ∈ [2.5, 3.5]

3. A_x sin(i)∈[60, 80] and P_orb ∈ [3, 4]

4. A_x sin(i)∈[80, 90] and P_orb ∈ [3, 4]

5. A_x sin(i)∈[90, 110] and P_orb ∈ [3.5, 4.0]

6. A_x sin(i)∈[110, 120] and P_orb ∈ [3.5, 4.0]

For each, three searches were performed in different frequency ranges: [0.333, 0.42], [0.366, 0.367], and [0.375, 0.385]. The |Ṗ/P| values range from 10⁻⁵ to 10⁻¹².

Each search stopped after 100,000 evaluations. However, if the current best power was found in the last 1,000 evaluations, the algorithm extended the search by a additional 10,000 evaluations. This process continued in increments of 10,000 evaluations until no improvements were detected in the last 1,000 iterations. As for the algorithm hyperparameters, the number of particles was set to n = 40, the cognitive component was set to c_p = 1.496172, the social component was set to c_g = 1.496172, and the inertia was set to ω = 0.72984, as in Bratton & Kennedy (2007).

Appendix B: additional table

All Tables

All Figures

Fig. 1. Light curve of X7 in the 0.3–10 keV energy band when using all available X-ray observations. The XMM-Newton, NuSTAR, Chandra, and Swift/XRT luminosities are not corrected for absorption and are taken as the average in each observation (assuming a distance of 7.5 Mpc). Chandra data are taken from Pintore et al. (2021). The inset shows the light curve after 2018 with the super-position of the profile of the tentative ∼240 d super-orbital periodicity (magenta).

In the text

	Fig. 2. Top: Lomb-Scargle periodograms of the Swift/XRT observations in the time intervals between 2020 and January 2023 in the 0.3–10 keV, 0.3–1 keV, and 1–10 keV energy range. Dashed lines indicate signals with > 3σ significance. Bottom: Full 1–10 keV Swift/XRT light curve folded with a period of 240 d.
In the text

Fig. 3. Background-subtracted light curve of X7 in XMM0, XMM1, and XMM2 (from left to right), binned at ΔT = 300 s (see the top-left legend), for the 0.3–10 keV (top panel), 0.3–0.8 keV (centre-top), and 0.8–10 keV (centre-bottom) energy bands. We also show the hardness ratio between the hard and the soft band (bottom panel). The dashed vertical lines divide the data between XMM0, XMM1, and XMM2, where the real time gaps between them were removed for displaying purposes only.

In the text

Fig. 4. Estimates of the X7 orbital period (Porb) and the projected semi-axis of the orbit for a spin period of ∼2.73s in XMM2. The confidence contour, more similar to an atoll rather than a ‘banana’ shape, is indicative of the low significance of the spin period detection. However, a similar plot can be obtained for XMM1 as well, boosting the real significance of the signal. The colour bar indicates the Leahy power of the signal.

In the text

	Fig. 5. Parameter space of the orbital parameters for the signal at 2.63s and 2.73s for XMM1 (left) and XMM2 (right), respectively. For displaying purposes only, in each plot the contours indicate the solution of the other observation. The colour bar indicates the Leahy power of the signal.
In the text

Fig. 6. Top: Hardness-intensity diagram for PN+MOS light curves of all XMM-Newton observations (after smoothing the data, see text). The data of the XMM0, XMM1 and XMM2 observations are reported in blue, orange and green, respectively. Hardness is defined as the ratio between the bands 2.0–10 keV (H) and 0.3–0.8 keV (S), while the rate is the sum of both bands. Smoothing with a moving window has been applied to the data (see text), and therefore we do not plot the errors on each point. The coloured and numbered segments on the top indicate the hardness ranges used for the spectra extraction shown in Fig. 8. Bottom: Hardness-intensity diagram for the available Swift/XRT data. Hardness is the ratio between the bands 1.0–10 (H) keV and 0.3–1.0 (S) keV.

In the text

	Fig. 7. Combined EPIC-pn + MOS(1,2) spectra of the six spectral states extracted from the HID (see text) unfolded through a power law of index zero. The colour convention is the same as indicated in Fig. 6-top.
In the text

Fig. 8. Combined EPIC-pn + MOS(1,2) spectra of the six spectral states extracted from the HID (see text and Fig. 7). The colour convention is the same as indicated in Fig. 6-top. The best fit is given by the TBABS(DISKBB + DISKPBB + CUTOFFPL) model (see Table 2), where the DISKBB, DISKPBB and CUTOFFPL models are indicated with dashed, dot-dashed and dotted lines, respectively. In the bottom panel of each plot, it is shown the ratio of data and model.

In the text

Fig. 9. Thermal component unabsorbed bolometric luminosities (0.001–30 keV) versus the soft (blue points) and hard (orange points) disc temperature (left), and thermal component temperatures (soft disc in blue, hard disc in orange) versus their inner disc radii (right) for Model 1. We assumed an inclination of the system of 10 degrees. The circle, triangle, diamond, cross, star, and square points indicate the spectra from one to six indicated in Table 1, respectively.

In the text

	Fig. 10. Thermal component unabsorbed bolometric luminosities (0.001–30 keV) versus the soft (blue points) and hard (orange points) disc temperature for Model 2. The circle, triangle, diamond, cross, star and square points indicate the spectra from one to six indicated in Table 1, respectively.
In the text

	Fig. 11. Combined RGS first and second order (grey and blue), EPIC-pn and MOS 1,2 (black, red, and green) spectra (top panel). The black and red lines refer to the continuum and continuum plus CIE models, respectively. The corresponding residuals are shown in the middle and bottom panels. We note that the Ne K – Fe L emission complex is near 1 keV.
In the text

	Fig. 12. Gaussian line scan for the combined RGS (0.4–2 keV) and EPIC (2–10 keV) spectra and different line widths. Rest-frame centroids of the typical dominant H-like and He-like transitions are labelled.
In the text

Fig. 13. Top panel: Multi-dimensional grids of emission models of plasma in collisionally ionisation equilibrium for the combined RGS (0.4–2 keV) and EPIC (2–10 keV) spectra. Middle and bottom panels: Grids of models of photoionised plasma in emission and absorption, respectively. Both scans adopted a line broadening of 100 km s−1 and solar abundances. The best-fit CIE model is also shown in Fig. 11.

In the text