NGC 2403 XMM4: evidence for a super-Eddington neutron star with a possible transient pulsation

We present a study of the X-ray source NGC 2403 XMM4 (4XMM J073702.2+653934) based on 20 years of archival observations with XMM-Newton, Chandra, Swift and NuSTAR. Although it has previously been classified as an ultraluminous X-ray source (ULX), we show that its luminosity rarely, if ever, passes the $10^{39} \rm ~erg~s^{-1}$ threshold luminosity for a ULX. It does, however, behave very similarly to ULXs, with its 0.3-10 keV spectra well described by two thermal components, the softer of which behaves consistently with the expectations for an advection-dominated disc ($L\propto T^{1.49\pm0.85}$), and we find tentative evidence for an extra spectral component above 10 keV. We also find moderately significant evidence for an absorption feature in one spectrum that could originate in an outflowing wind, although a cyclotron resonance scattering feature is also a possibility. Most intriguingly, we find a possible transient pulsation at $\sim 3.32$ Hz in a short segment of one observation using an accelerated pulsation search. This evidence suggests that NGC 2403 XMM4 is displaying many of the hallmarks of super-Eddington accretion at luminosities between $5 - 10 \times 10^{38} \rm ~erg~s^{-1}$ which, when considered alongside the putative pulsation, points to the presence of a neutron star as the accreting object this system.


INTRODUCTION
Ultraluminous X-ray sources (ULXs) are defined as point-like X-ray sources with typical X-ray luminosities ( X ) > 10 39 erg s −1 that are located in the extra-nuclear regions of galaxies.Historically, the two main interpretations for such a relatively high luminosity have been either accretion at sub-Eddington rates onto a massive compact object such as an intermediate mass black hole (IMBH), or super-Eddington (or ∼Eddington) accretion onto a less massive object, e.g. a stellar mass black hole (BH) or even a neutron star (NS).ULXs are now becoming increasingly well studied, having been subject to dedicated observations for more than two decades (see e.g.Kaaret et al. 2017;King et al. 2023 for recent reviews).One surprising result in the field was the first detection of a pulsating signal, from the ULX M82 X-2 (Bachetti et al. 2014) which provided the first strong evidence that ULXs could be powered by super-Eddington NSs, i.e.NS-ULXs, in sharp contrast to the original assumption that they would most likely be powered by BHs.Following this first detection of a pulsating NS-ULX -hereafter PULX -several more objects were revealed to have a pulsating signal in their observational data, with notable examples including NGC 5907 ULX-1 (Israel et al. 2017a), NGC 7793 P13 (Fürst et al. 2016;Israel et al. 2017b), M51 ULX-7 (Rodríguez Castillo et al. 2020) and NGC 1313 X-2 (Sathyaprakash et al. 2019); see also Table 2 of King et al. (2023) ★ E-mail: wasutep@g.swu.ac.th (WL) and references therein for other detected PULXs.In addition, NS-ULXs can also potentially be identified by the detection of cyclotron resonance scattering features (CRSF) imprinted on ULX spectra, due to the presence of strong magnetic fields (≳ 10 12 G) in these objects (e.g., Brightman et al. 2018;Walton et al. 2018b).Indeed, there is now the possibility that a majority (or, at least, a substantial fraction) of the whole ULX population could be powered by NSs, (c.f.Walton et al. 2018a) which implies by the observational definition of ULX luminosity that these NSs must be highly super-Eddington accretors.
NGC 2403 is a nearby, Scd spiral galaxy located at distance of 4.2 Mpc (Tully & Fisher 1988), that provides a good opportunity for studying the brightest X-ray sources residing in a normal galaxy.Indeed, the X-ray sources in this galaxy have been well studied with nearly 200 discrete sources having been identified (Schlegel & Pannuti 2003;Binder et al. 2015).Among these, there is one source -CXOU J073625.5+653540-that has been shown to be a persistent ULX and has been analysed by several studies using data from various observatories, e.g.Kotoku et al. (2000); Liu & Bregman (2005); Winter et al. (2006); Isobe et al. (2009).Furthermore, the detection of a transient X-ray pulsar (CXOU J073709.1+653544) in this galaxy, with a maximum luminosity of 2.6 ×10 38 erg s −1 , has also been reported by Trudolyubov et al. (2007) in which it was argued that this luminosity exceeds the isotropic Eddington limit for a 1.4M ⊙ compact object.Given that, there are at least two strong candidates for super-Eddington sources identified to date in this galaxy.
In fact, there is another X-ray source that has been classified as a ULX candidate in NGC 2403, in the catalogue published by Walton et al. (2011): 2XMM J073702.1+653935.1This is a relatively low luminosity ULX candidate, in the poorly-studied 'Eddington threshold' regime discussed in Earnshaw & Roberts (2017) ( X < ∼ 10 39 erg s −1 ).It was previously studied in several papers, where it was referred to by a variety of different names: e.g., NGC 2403 X3 (Liu & Bregman 2005), NGC 2403 XMM4 (Winter et al. 2006), NGC 2403 source 5 (Kotoku et al. 2000;Isobe et al. 2009) etc. (hereafter, we will refer to this source as NGC 2403 XMM4 for convenience).In fact, most works that studied NGC 2403 XMM4 were actually population studies and so the source was analysed together with, and then compared to, numerous other sources (see e.g.Schlegel & Pannuti 2003;Liu & Bregman 2005;Winter et al. 2006).These previous studies found similar results: that the source could be well approximated by an absorbed, single component, hard power-law model (with photon index Γ ∼1.6 -1.9 ) or an absorbed thermal bremsstrahlung model with the plasma temperature of ∼5 keV, while its luminosity was reported to be ∼3 -5 ×10 38 erg s −1 .Given the source spectral shape and luminosity, it was proposed by Winter et al. (2006) that NGC 2403 XMM4 could be classified as a sub-Eddington intermediate mass BH.However, in-depth analysis of this source is still required to firmly identify the class of compact object powering it.
Given that NGC 2403 has been observed multiple times by various observatories since the beginning of this century, including several observations in the current decade, there is now sufficient archival data to provide us with the opportunity to analyse this source in detail.In this work, we provide a comprehensive analysis of the Xray spectral and temporal characteristics of NGC 2403 XMM4, with the aim of investigating whether can we confirm that it is another super-Eddington source in NGC 2403, and if so whether we can determine the nature of its compact object.The paper is laid out as follows.In Section 2, we detail the observational data used in this study and explain how we produce the science products for analysis.The analysis methods and results are shown in Section 3 and these are then discussed in Section 4, and our conclusions made in Section 5.

OBSERVATIONS AND DATA REDUCTION
In order to study the properties of NGC 2403 XMM4, we gather observational X-ray data from the archives of the current generation of X-ray observatories -i.e.XMM-Newton, NuSTAR, Chandra and Swift -that have sufficient data quality for timing and spectra analyses.The data used in this work are shown in Table 1, excepting the Swift data which are available in Appendix A. The detail of observations analysed in this work and how they were reprocessed and converted into science products are described below.

XMM-Newton observations
We collected all XMM-Newton archival observations that contain data from the X-ray source NGC 2403 XMM4; the data were downloaded from the XMM-Newton Science Archive (XSA).2However, the XMM-Newton observations analysed in this work are ones in which the summation of good exposure times (defined as the good time intervals, or GTIs) for each detector after the removal of background flaring periods is > ∼ 10 ks, to ensure that we obtained sufficient data quality for our analysis.Using this criterion, we ended up Note. The total (summed) value of good time intervals (GTI) for each instrument; for XMM-Newton, it is the total GTI of pn, MOS1 and MOS2 detectors, respectively, while for NuSTAR, it is that of FPMA and FPMB detectors, respectively. The background-subtracted, instrument count rate of the source; the value is obtained from pn and ACIS-S detectors in the 0.3 -10 keV band for XMM-Newton and Chandra, respectively, while for NuSTAR, the 3.0-30 keV count rate obtained from the FPMA detector is reported.† The number in the square bracket indicates the epoch of data used in the broadband spectral analysis of XMM-Newton + NuSTAR data (see Section 3.1.4).
with the five XMM-Newton observations presented in Table 1.However, we note that there are another five observations in which the total GTIs are insufficient for analysing in depth (i.e., Obs ID 0150651101, 0150651201, 0301860301, 0729560901 and 0870870301); here we were, at least, able to measure and report their individual fluxes, as shown in Fig. 1, using a similar method to that we used to measure the flux of individual Swift observations (see Section 3.1).The observational data were reprocessed following the standard procedure suggested on the XMM-Newton data analysis webpage 3 using the Science Analysis System (sas) software version 21.0.0 with the latest calibration files at the time that the data were reprocessed (May 2023).In detail, we reprocessed the pn and MOS data using the SAS packages epproc and emproc, respectively.The reprocessed files were then inspected for any periods that were affected by background flaring events by visually checking the high energy light curves in the 10-12 keV energy band.If there were any observing periods in which the count rates exceeded > ∼ 0.4 count s −1 for the pn detector or > ∼ 0.35 count s −1 for MOS1 or MOS2, then the data from these periods were removed from the reprocessed event files.For spectral analysis, we created the spectral files and their corresponding redistribution matrix file (RMFs) and auxiliary response files (ARFs) for all observations using the default event flags and patterns for each detector, using a circular data extraction region with radius 25 ′′ located at the source position.We used a source-free region of 32 ′′ radius located nearby and on the same detector chip as the source to extract a background spectrum.We note that for the XMM-Newton spectra that we analysed simultaneously with NuSTAR data, we also applied a correction to the effective areas in the ARFs by setting the parameter applyabsfluxcorr in the SAS task arfgen to 'yes'.This could help to avoid a flux discrepancy between XMM-Newton data and external flux data (i.e., those of NuSTAR). 4inally, all spectra were grouped by the SAS task specgroup setting the parameters mincounts = 20 and oversample = 3, to utilise the chi-squared statistic in spectral modelling.
We also prepared files for timing analysis, with the goal of searching for pulsations.We first performed a barycentric correction to the cleaned event files using the SAS command barycen, and then created both pn and MOS light curves for all observations using the same source and background regions as used in the spectral extraction, and all science events with the appropriate flags.The binning of the resulting light curves was set to be equal to the temporal resolution of the science mode of each instrument used during the observation, i.e., 73.4 ms for the pn and 2.6 s for MOS1 and MOS2.We also extracted the light curves in three separate energy bands: 0.3 -2 keV (soft band), 2 -10 keV (hard band) and 0.3 -10 keV (full band).

NuSTAR observations
The three NuSTAR archival observations were reprocessed using the task nupipeline of NuSTAR Data Analysis Software (nustardas) version 2.1.2,which is part of the heasoft package. 5We also applied a barycentric correction to the cleaned event files using the ftool task barycorr. 6Then, spectra and timing products were extracted from the cleaned event files using the task nuproducts; the source region for extracting science products was defined as a circular region with a radius of 40 ′′ located at the source position, while the background was taken from a circular, source-free region with a radius of ∼120 ′′ -150 ′′ close to the source region.Similar to XMM-Newton spectra, the NuSTAR spectra were then grouped to have a minimum of 20 counts per bin for spectral analysis.

Chandra observations
There are four Chandra observations in which NGC 2403 XMM4 lies inside the field of view, three of which were taken in 2004 and are ToO observations of supernova 2004dj in NGC 2403 (Pooley & Lewin 2004).All Chandra data were download from the Chandra data archive7 and then reprocessed using the standard script chan-dra_repro which is part of the Chandra Interactive Analysis of Observations (ciao) software version 4.15.8Data were checked for background flaring, which was not significant in 3/4 observations, and only minimally affected the fourth (Obs.ID 2014, from which < ∼ 2 per cent of data was removed using the task lc_sigma_clip.In addition, we also applied a barycentric correction to the time column in the cleaned event files using the ciao tool axbary. We then extracted science products for all observations in the full (0.3-10 keV) band using the ciao task srcflux; this task produced both spectral and timing products (i.e., light curves) and provided these for each individual dataset as well as for merged data.Here, the source data were extracted from a circle of 5 ′′ radius, while the background region was a source-free annulus of 7 ′′ inner radius and 12 ′′ outer radius; both regions were centred on the source location.Note that for the light curves, the time resolution that we used for binning the data is set equal to the nominal frame time of the Chandra ACIS detectors in the observations used in this work, i.e. ∼3.2 seconds.

Swift observations
The Swift data for NGC 2403 XMM4 were selected from observations in the Swift data archive in which the ULX is within the field of view.In total, we obtained 27 useful observations that have sufficient net photon counts to constrain the ULX count rate and flux; the observational details are summarised and shown in Table A1 in Appendix A.
The Swift data were reprocessed to create level 2 cleaned eventlist files using the pipeline xrtpipeline version 0.13.7, which is included in the heasoft software package.Then, the spectrum of each observation was extracted from a circular aperture with radius of 35 ′′ centred on the ULX position.A background spectrum was obtained from a source-free annulus around the source position, which has 50 ′′ inner and 100 ′′ outer radii.The detail of the Swift spectral analysis will be described in Section 3.1.

DATA ANALYSIS AND RESULTS
In this section we present the results from spectral and temporal analyses of the data from NGC 2043 XMM4.As our investigation is largely focussed on whether this object is a super-Eddington accretor, given its previous identification as a ULX candidate, we will compare its properties to those of ULXs throughout this section.

Spectral analysis
Throughout this paper, spectral analysis was performed using the Xray spectral fitting package xspec, version 12.13.0(Arnaud 1996). 9e also assume that the observed X-ray data are absorbed by material along the line of sight in the direction to NGC 2403 XMM4 so that two absorption components (tbabs×tbabs in xspec) were added to every spectral model we used in this paper.The first component represents the absorption by neutral material in our own Galaxy, for which the value was frozen at 0.04 × 10 22 atoms cm −2 (Dickey & Lockman 1990),10 , while the other component can be interpreted as the absorption by material in the host galaxy and/or by material in the close proximity of the source itself, and we refer to this as the external absorption column density hereafter.If not otherwise specified, we allowed the value of the latter component to be a free parameter during the spectral fitting, and the spectra were normally fitted in a range of 0.3 -10 keV.The source luminosities were calculated assuming a distance of 4.2 Mpc (Tully & Fisher 1988).In addition, for those observatories where multiple detectors take simultaneous data, a constant component was added to the models used to fit the data to allow for any small differences between detector responses.The errors on spectral fits are quoted as 90 per cent confidence ranges of parameters unless otherwise specified.

Single component models
We began the analysis by fitting the spectra of individual Swift observations.However, due to the low number of net photons of the source (∼5 -50 counts per observation), the xspec statistic cstat was applied to fit the Swift spectra, that was kept unbinned for this analysis.We also simplified the spectral fitting by using a single component absorbed power-law model and reduced the number of free parameters by freezing the absorption column density and power-law photon index at 0.2 × 10 22 atoms cm −2 and 1.6, respectively; note that the fixed values of these two parameters were obtained from the best-fit values of the merged Swift spectrum fitted by the same model (see Table 2).Moreover, these parameter values were also consistent with those we could obtain from individual XMM-Newton spectra fitted by an absorbed power-law model; this implies that the parameter values are a good (albeit rough) representation of the power-law-like spectral shape of this source.Given that the two parameters are frozen, the only free parameter in the fitting is the power-law normalisation which was calculated in terms of flux using the xspec cflux convolution component and then converted into luminosity.Fig. 1 shows the long term luminosity of NGC2403 XMM4 as observed by Swift in the interval 2006 -2021, and augmented by single epoch flux measurements from the other detectors considered.
To investigate the Swift spectra further, we maximised the spectral quality by merging all the individual spectra to become a single averaged spectrum using the ftools task addspec11 and grouped the data to have a minimum of 20 counts per bin.We then fitted the spectrum with two traditional, single component models widely used for ULXs, an absorbed power-law and an absorbed multi-colour disc blackbody (MCD; diskbb in xspec).The result is shown in Table 2 in which both models are seen to fit to the data equally well, while the best-fitting values of the MCD temperature and power-law photon index are consistent with past studies of ULXs (e.g.Swartz et al. 2004;Gladstone et al. 2009).
Another point we would emphasise is that the averaged luminosity of NGC 2403 XMM4 obtained from Swift data is ∼0.6 -0.8 × 10 39 erg s −1 , which lies below the traditional threshold of ULX luminosity, although individual data points show that the source occasionally rises to luminosities above 10 39 erg s −1 .Given this behaviour, we discuss whether the source should be regarded as a ULX in Section 4.1.
Next we analysed the Chandra data.Given these observations were taken in the period 2001-2004 we note this extends coverage of NGC 2043 XMM4 to earlier times than the Swift data.We began by analysing each individual Chandra dataset in the 0.5 -7 keV energy band, with the spectra binned to have a minimum of 15 counts per bin, and used the same absorbed, single component models as per the Swift merged data.We found that although both models could yield acceptable fitting results, the absorbed MCD model is able to describe the individual Chandra spectra better (a total bestfitting  2 of 317 for the MCD model versus 344 for the powerlaw, for a total of 309 degrees of freedom).The calculated fluxes from the MCD model are shown in Fig. 1 which indicates that the flux obtained from Chandra data is roughly constant.Furthermore, the best-fitting MCD temperatures of all the spectra are consistent, with their values ≈1.9±0.2 keV, which could imply that they are all in a similar spectral state.Interestingly, the temperatures for the disc component are relatively high when compared to stellar-mass Galactic BHs in the thermal dominant state (c.f.fig. 1 of Miller et al. 2004; see also Done et al. 2007 and references therein), however similar high disc temperatures can sometimes be seen for Z-type NS low-mass X-ray binaries (e.g., Mondal et al. 2018;Homan et al. 2018).Similar high disc temperatures are also seen in the broadened disc-like regime (Sutton et al. 2013, see also Luangtip et al. 2016), that is sometimes observed for ULXs.
Given the obvious similarity of the individual Chandra spectra, we decided to improve the spectral data quality by merging all the spectra together.The improvement in signal to noise was significant enough that the data was appropriate for fitting with more complicated double component models, which we detail in the next sub-section.

Double component model
The individual XMM-Newton spectra and merged Chandra spectrum have sufficient data quality to enable the fitting of more complex models than the single component fits in the previous section.Here, we attempt two component fits based on two thermal components, as is now widely accepted in ULX spectroscopy (Stobbart et al. 2006;Koliopanos et al. 2017;Walton et al. 2020), although the physical interpretation of these components varies.We choose to model the data using a blackbody continuum and a MCD -tbabs×tbabs(bbody + diskbb) in xspec -where, for example, the soft component may be the optically thick base of an outflowing wind (as detected by e.g.Pinto et al. 2016), and the harder component may come from the innermost regions of the accretion flow.The fitting results are shown in Fig. 2 and best-fitting values are shown in the top panel of Table 3; the calculated best-fit luminosities are also displayed in Fig. 1 to further enhance the long timescale light curve of NGC 2403 XMM4.The model provides an excellent fit to the data, with reduced  2 ( 2  =  2 /dof) values below unity in 5 out of 6 cases.However, the spectrum of XMM-Newton observation 0870870501 has a slightly higher, but statistically acceptable, value of  2  = 1.13, corresponding to a p-value for the null hypothesis probability of 0.2.In fact, the spectral data from this observation appear to exhibit structured fitting residuals at ∼1 keV and ∼1.7 keV (see Fig. 2 middle-right image); we will return to this point in section 3.1.3.
In addition, recent work on several ULXs (e.g.Robba et al. 2021;Barra et al. 2022) has shown correlations between the component fluxes and temperatures in similar two-component models.In fact, starting from Miller et al. (2013) and the following work (e.g.Middleton et al. 2015a;Walton et al. 2017), the correlations, especially in the case of the soft component, tend to be detected when ULXs do not show strong variability in their absorption column.Here, we also check if this is the case for this source.Excepting for the dataset 0870870201, the best-fit column density values are consistent within errors.So, excluding that observation, we checked if any correlation is present by repeating our spectral fits with the column density held fixed across the five remaining datasets.A global  H of 0.05 ± 0.02 × 10 22 atom cm −2 was obtained (see lower part of Table 3), while the best-fit spectra and the relation between the temperature and the intrinsic flux of both soft and hard components are plotted in Fig. 3. Indeed, some degree of positive correlation between the soft component temperature and the flux is suggested in the plots although a Pearson's correlation () measurement shows that the correlation is not statistically significant.However, this could be driven by a second outlier dataset, 0870870801, with a significantly lower flux in its hard component, and a best-fit spectral shape suggestive of two distinct components, similar in form to the hard ultraluminous regime as defined by Sutton et al. (2013).Interestingly, if we also remove the outlier dataset from the correlation measurement, a linear correlation for the remaining four datasets is detected using a Pearson  correlation ( = 0.94), with a best-fit slope of 1.49 ± 0.85 (1 error), implying a luminosity () -temperature () correlation for the BB component in which  ∝  1.49±0.85 .

Narrow absorption-like residuals
As we mentioned above, the spectra appear to have some structured residuals to the best-fitting continuum model in XMM-Newton observation 0870870501, most notably at ∼1 keV and ∼1.7 keV.We investigated this further by adding multiplicative gaussian absorption line components (gabs in xspec) to the model.Interesting, the additional components do help to improve the fitting; the left panel of Fig. 4 shows the new best-fitting model together with the fitting residuals in which the best-fitting absorption line energies are at 0.99 +0.05 −0.04 keV and 1.68 ± 0.06 keV.The addition of both components improves the fit, with the absorption component at 0.99 keV improving the fitting statistic by Δ 2 = 5.63,12 while that of the component at 1.68 keV changes the fitting statistic by Δ 2 = 17.38, for three fewer degrees of freedom in each case.Evidently, the putative absorption feature at 1.68 keV is relatively more significant than that at 0.99 keV.We also note that these residuals cannot be explained by either the complex, unresolved emission lines of a collisionally-ionised diffuse gas (apec model in xspec) or a warm absorber, as calculated by the XSTAR model (Kallman & Bautista 2001). 13ur next step was to attempt to place better constraints on the significance of the putative absorption features.We used the xspec command fakeit to simulate pn, MOS1 and MOS2 spectra using the same observational response files as used for the real data, the best-fitting two component model from the 0870870501 data shown in Table 3 (i.e., the two component model without the additional absorption line features), and an exposure time for each instrument    .49, while the grey dashed and grey dashed-dot lines outline the  ∝  2 and  ∝  4 correlations, respectively.Note that the plotted data are obtained from fitting with the  H tied, excepting for that of XMM-Newton 0870870201 in which the plotted data are adopted from the individual fit (see text for detail.) exactly matching the real data (see Table 1).In addition, realistic statistical fluctuations were also added to the simulated spectra by fakeit.In total 60,000 sets of simulated spectra were created.Each set of fake observational spectra was then analysed by the same models used to analyse the real data, i.e. (i) two thermal components and (ii) two thermal components with an additional absorption component model; note that for the latter model, the parameters of the absorption component, including the line energy, were allowed to be free parameters during the fitting, similarly to the method that was used to analyse the real data, and the improvements in  2 from the additional component were recorded.Hence, we are able to determine the probability of detecting features that are simply statistical fluctuations in the simulated data, rather than real characteristics of the source.We found that the probability of obtaining values of Δ 2 > 5.63 and >17.38 from statistical fluctuations alone are 16.28 per cent and 0.53 per cent respectively, corresponding to significance levels of ∼ 1.1 and ∼ 2.6.Hence, there is no strong evidence that either feature is real, although the feature at 1.68 keV is at least marginally significant.We will discuss the believability of this putative absorption feature and its possible physical origin further in Section 4.

Broad band spectral model
Near-simultaneous observations of the galaxy NGC 2403 were obtained with XMM-Newton and NuSTAR near to the end of October 2020 and April 2021.NGC 2403 XMM4 lies in the field of view of both sets of detectors, so that they provide us with the opportunity to study the broadband spectra of this source.However, there is an offset of a few days between the start of observations by the two telescopes in both epochs.Given this, we grouped together two sets of data that had an offset of ∼1-2 days, which we hereafter refer as epoch 1 and epoch 2 data (see Table 1).For the spectra in epoch 2 that have more than one observation per observatory, we combined all spectra and response files obtained from the same detector, using the sas task epicspeccombine and the ftools task addspec for XMM-Newton and NuSTAR data, respectively.Throughout the analysis in this section, XMM-Newton and NuSTAR spectra from each epoch were modelled simultaneously within the full energy range of 0.3-30 keV.
We began by fitting an absorbed, two thermal component model -exactly the same model as was used in Section 3.1.2-to examine whether a two component model is sufficient to explain the broadband spectra of the source, especially at the high energy end of the spectra.The best-fitting result is shown in Table 4, along with the spectra and residual plots in Fig. 5.The reduced  2 values in both epochs are ∼1, implying that the model is a reasonable description of the data.However, given both previous results for ULXs that show a hard excess above 10 keV (Walton et al. 2018a), and that, in particular, epoch 1 (Fig. 5 middle panel) shows positive residuals at ∼10 -20 keV in which the model seems to underestimate the data, this convinced us to investigate further by adding a third component to the model -an exponential cutoff power-law (cutoffpl in xspec)in an attempt to explain the spectra at high energies.The addition of this component led to a Δ 2 of ∼10 for epoch 1, and ∼3 for epoch 2, for three fewer degrees of freedom.This suggests that the additional component marginally improves the fit to the epoch 1 spectra, while there is little improvement in that of the epoch 2 spectra.Thus, to examine whether the cut-off power-law component helps to improve the fit significantly, the xspec statistical tool ftest, i.e., the F-statistic probability of the fit improving given the Δ 2 and the change in degrees of freedom before and after adding the extra component, was applied as an first order approximation14 .Note that the xspec ftest always calculates the probability in term of its p-value; here, the values of 1 -p are reported to better clarify the meaning of the probability.The result is shown in column 8 of Table 4.This does suggest that the epoch 1 spectra are likely to require the additional exponential power-law component to explain their high energy tail (with a probability of ∼98 per cent, corresponding to a confidence level of ∼2.3) while this component appears not to be needed for the epoch 2 spectra.We will discuss this point further in the discussion section.

Timing analysis
We present the long-term light curve of NGC 2043 XMM4 in Fig 1.This demonstrates that it is a remarkably persistent source on timescales of days to years, with all observations over a ∼ 20 year period having observed fluxes in the range ∼ 0.6−1.3×10 39erg s −1 .Here, we investigate whether this persistence is a characteristic down to short timescales, by analysing the timing information from the individual observations in Fourier space.We start by searching for any significant variability in the observed frequency domain using the ftool task powspec, using the extracted light curves in the energy bands for each instrument defined in Section 2, and applying the Leahy normalisation (Leahy et al. 1983) in which the Poisson noise power is set equal to two.Fig. 6 illustrates the power-spectral densities (PSDs) of XMM-Newton observation 0164560901 -which is one of our highest quality timing datasets -in the full, soft and hard energy bands.Unfortunately, there is no variability detected significantly above Poisson noise in any individual observation.Although the analysis was also applied to analyse the concatenated light curves of all XMM-Newton and Chandra data, similarly, no variability was detected in these light curves either.We double-checked these results using an alternative software package, namely hendrics version 7.0 (Bachetti & Huppenkothen 2018) 15 which is a set of command-line scripts based on the stingray Python library (Huppenkothen et al. 2019a,b;Bachetti et al. 2023).No significant variability was found in any PSD created by this task either.
Although there is no periodical modulation hinted at by the calculated PSDs, we cannot immediately rule out NGC 2403 XMM4 as an X-ray pulsar.PULXs do not often show strong pulsations in their PSDs as they tend to display strong period derivatives (  ), as a result primarily of orbital motion, that result in blurring of the pulsation power across multiple frequency bins (see e.g.Sathyaprakash et al. 2019).Thus, we began searching for a pulsating signal in the frequency range of 0.1 -5 Hz using data from the only instruments with sufficient temporal resolution to investigate this regime, the XMM-Newton pn and NuSTAR FPMA and FPMB detectors.We used  Notes. Photon index and  cut-off energy of the exponential cut-off power-law component. The probability indicating that the data would require the third cutoff power-law component to explain the spectra (see text for detail).* No error bars reported as the parameter was fixed at its best-fit value to enable constraints to be placed on the confidence intervals of the other free parameters.
The rest of the column definitions are as per Table 3.
the full energy band of the data (i.e., 0.1 -12 keV for pn data and 3 -79 keV for FPMA & FPMB) and analysed it using the hendrics task henaccelsearch which utilises the accelerated search method detailed in Ransom et al. (2002).Subsequently, any detected candidates were further analysed and verified by the hendrics epoch folding algorithm henzsearch which applies the  2 statistic with  harmonics (c.f.Buccheri et al. 1983); here =3 (i.e.,  2 3 ) was used.Using the above method, the search was applied to the whole length of each individual observation, but no significant pulsation could be verified in any searches.However, as illustrated in e.g.Sathyaprakash et al. (2019), the source might not exhibit a persistent observable pulsating signal throughout the entire observational period and, indeed, the signal might not be detected in all observational data; instead, the pulsating power might only be strong enough to be detected during some discrete periods.To check this, we divided each observational dataset into a number of segments.Instead of dividing each observation into the segments with equal time length, we used the gap between GTIs -i.e., observational periods with bad background flaring or no data -as boundaries for segmenting each dataset.This is to ensure that the gaps in the data are excluded from the calculations.Fig. 7 illustrates how the observational data were divided into segments using this method; note that, while most segments are sufficiently long to be useful for the calculations, some segments are too short (∼0.1-0.3 ks) and so were excluded from the analysis.
By dividing the observational data into segments and applying the search again, fascinatingly, there is one time segment in the XMM-Newton data of observation 0164560901 (observed in September 2004) which exhibited a detectable pulsating signal -i.e. the fourth segment in Fig. 7 -at ∼ 3.32 Hz (i.e., ∼0.3 s period) with a  2 rather discrete; it is not obviously sinusoidal in profile, as per most PULXs (e.g.Bachetti et al. 2014).In addition, we also plot the  - parameter space in Figure 9 to demonstrate the peak location where the statistical significance of these two parameters are highest; the analysis constrains values of  of 3.3188 +1.95×10 −6 −2.26×10 −6 Hz and the  of -1.61 +0.0574 −0.0240 ×10 −7 Hz s −1 .Note that the range over which the dataset segment was searched was  = 3 -3.5 Hz with resolution of ∼ 3.95 × 10 −5 Hz, while that of  was ∼(−2.46 -+2.46) × 10 −7 Hz s −1 with resolution of ∼2.14 × 10 −8 (i.e., 291,456 trials in total).We emphasise that there is no other significant pulsation detected in any other segment of this observation, or in any segments of the other observations, including those of the NuSTAR data.However, we do note that there were at least two lower significance detections with  2 3 powers of < ∼ 40, detected at frequencies of ∼3.2 Hz and ∼3.4 Hz in segments of the data from XMM-Newton observations 0870870501 and 0870870801.Although it is very encouraging that we detect further signals in close proximity to the frequency of our one significant detection (see Sathyaprakash et al. 2019 for the importance of such detections in confirming real, weak pulsations), their low statistical significance is insufficient in either case to confirm they constitute real pulsating signals.

DISCUSSION
In this work, we have studied observational data for the X-ray source NGC 2403 XMM4 obtained from the archives of multiple X-ray telescopes from the current generation of observatories.Combining the variety of data we obtained provides a broadband energy range from soft to hard X-rays (i.e., ∼ 0.3 -30 keV), and temporal resolution ranging from as short as ∼ 5 Hz to long timescales of up to two decades.This range of data permits us to conduct the variety of spectral and temporal analysis shown in the previous section.Here, we discuss what they mean for the nature of this object.

Should we regard this source as a ULX?
Our interest in NGC 2403 XMM4 originates from its inclusion in several catalogues and/or studies of ULXs (see e.g.Walton et al. 2011;Winter et al. 2006).However, given the 20 year light curve of the source we show in Fig. 1, it is obvious that, most of the time, the source luminosity lies below 10 39 erg s −1 , and hence is below the threshold luminosity normally used to define ULXs.We therefore begin by posing a question: should we continue to regard this object as a ULX?
The only data points in Fig. 1 consistent with a luminosity above the threshold for ULXs were obtained from Swift observations, and all of those may be consistent with luminosities below this, given the low number of counts per Swift observation and hence relatively high uncertainties.Indeed, we cannot reject the null hypothesis that the source luminosity is persistently below 10 39 erg s −1 at 3 confidence from this dataset.Furthermore, there is a possibility that we have overestimated the distance to NGC 2403.Our value of 4.2 Mpc -which is adopted from the Walton et al. (2011) catalogue -is actually at the higher end of the range of measurements for the distance to NGC 2403, as reported in the NASA/IPAC Extragalactic Database (NED) 16 , which typically range from ∼ 3 − 4 Mpc (e.g. de Vaucouleurs et al. 1981;Bottinelli et al. 1984;Olivares et al. 2010).Thus, if the source distance is reduced by a factor of 3/4, the luminosity would be reduced by a factor of (3/4) 2 ∼ 0.6.This would imply that the source luminosity lies well below the traditional ULX luminosity regime, at ∼ 5 × 10 38 erg s −1 .Indeed, this is very likely why this source is included in the Walton et al. (2011) catalogue, but not in more recent ULX catalogues (e.g.Earnshaw et al. 2019;Walton et al. 2022) despite being detected in multiple serendipitous surveys.
It might however be argued that 'ultraluminous' is a better description of a behaviour of an object, than of the nature of the object itself; indeed it has long been suggested that ULXs represent a particular behaviour of stellar-mass compact objects (e.g. the 'ultraluminous state' of Gladstone et al. 2009).There is now evidence of many objects catalogued as ULXs occasionally (or in some cases frequently) dropping to luminosities below the ULX regime, for example some of the Eddington threshold objects studied by Earnshaw & Roberts (2017), or the sample of highly variable ULXs presented by Song et al. (2020) that notably highlights that several of the known PULXs do this.Hence, it is still plausible that NGC 2403 XMM4 is a ULX, at least some of the time, although the evidence based solely on its luminosity is somewhat tenuous.

Does NGC 2403 XMM4 behave like a ULX?
The NGC 2403 XMM4 spectra can be modelled by two-component models composed of two thermal components, consistent with ULXs (e.g.Stobbart et al. 2006;Koliopanos et al. 2017), as is illustrated by Fig. 2 and Fig. 3.The spectra of this source appear disc-like at higher luminosities, but display two distinct components at lower luminosity.By comparison to the spectral regimes proposed by Sutton et al. (2013) for ULXs, the former and latter spectral shapes can be classified as broadened disc and hard ultraluminous, respectively.In their study, Sutton et al. (2013) showed that the dominant spectra of ULXs (as a population) changed from broadened disc below ∼ 3 × 10 39 erg s −1 to hard (or soft) ultraluminous as the luminosity increases, which they suggested might be an evolution from sources accreting at around the Eddington rate to overtly super-Eddington accretion flows.However, subsequent work demonstrated that the spectra of BHBs in the thermal-dominant state, typically accreting matter at ∼10 -30 per cent of their Eddington rate (and so having luminosity ∼ 10 38 erg s −1 ), exhibit a broadened-disc like spectral shape, similar to that of the low-luminosity ULX spectra, in the 0.3-10 keV band that ULXs are observed in (Sutton et al. 2017).Hence, this spectrum could suggest that the source is a BH binary accreting at moderately high rates, but still in the sub-Eddington regime.However, studies of individual ULXs (see e.g.Luangtip et al. 2016) have demonstrated that these objects evolve from a distinct twocomponent (hard ultraluminous) spectrum at lower luminosities, to a single broad disc-like component (broadened disc) at the highest luminosities.We cannot rule out that that is happening for NGC 2403 XMM4, despite it much lower luminosity.
A further hint that NGC 2403 XMM4 is accreting at super-Eddington rates comes from broadband spectroscopy.The presence of a steep high energy tail, albeit only at moderate statistical significance (∼ 2.3) in one of two observational epochs, is another similarity to ULX spectra.Such a component is a general property of ULXs and PULXs, in addition to the two thermal components below 10 keV (e.g.Mukherjee et al. 2015;Walton et al. 2017;Koliopanos et al. 2019;Walton et al. 2020).This high energy tail can be explained in the context of Compton up-scattering in an X-ray corona (a similar mechanism to that found in sub-Eddington BH binaries) or a super-Eddington accretion column onto a magnetised NS (Walton et al. 2020); the former scenario could be possible for both BH and non-magnetic NS accretors whilst the latter requires that the accretor must be a magnetised NS.
Moreover, after careful removal of the outlier datasets, we found a correlation between the flux and temperature of the soft component of  ∝  1.49±0.85 .This result is broadly consistent with those found in previous ULX studies, including for PULXs (e.g.Pintore & Zampieri 2012;Walton et al. 2013Walton et al. , 2020;;Robba et al. 2021;Barra et al. 2022).Although many of these these studies detect a putative  − correlation for both soft and hard components, in our case a correlation is only found for the soft component, which may be related to the small number of data points and their relatively low spread of values in the warm component.Interestingly, the temperature power index of ∼1.5 is consistent with one report of a similar relationship for the PULX NGC 1313 X-2, in which an index of 1.2 ± 0.3 was reported by Pintore & Zampieri (2012) despite leaving  H free to vary; however a subsequent analysis of a larger dataset for this source found a very different  ∝  −3.9±1.0 relationship (Barra et al. 2022).The relationship we derive for NGC 2403 XMM4 is consistent within errors with the  ∝  2 prediction for advection-dominated (slim) disc models (Watarai et al. 2000;Watarai & Mineshige 2001), that may form at the centre of super-Eddington accretion flows (Poutanen et al. 2007).It is only very marginally consistent with a standard disc model where  ∝  4 (Shakura & Sunyaev 1973), at the ∼ 2.95 level.This behaviour therefore also points to the possibility that the accretion in NGC 2043 XMM4 is super-Eddington.

Is this a super-Eddington neutron star?
Although NGC 2403 XMM4 appears to spend little time (if any) above the 10 39 erg s −1 threshold for a ULX, its spectral behaviour does appear distinctly ULX-like, which may imply that it is accreting at super-Eddington rates.If so, we may estimate its compact object mass simply by using equation 1 of Feng & Soria (2011; see also Frank et al. 2002), assuming that the time averaged luminosity of NGC 2403 XMM4 is ∼5 × 10 38 erg s −1 and the emission is isotropic.This luminosity then corresponds to a compact object mass <3.8M ⊙ for an object accreting matter at super-Eddington rates, or even lower if the X-ray emission is beamed.This implies that the compact object is a NS; for a theoretical mass range for NSs of 1.2 -2.9M ⊙ (Kalogera & Baym 1996;Suwa et al. 2018), the typical luminosity of ∼5 × 10 38 erg s −1 corresponds to accretion at between 1.3 -3.2 times Eddington; for peak luminosities up to ∼ 10 39 erg s −1 , this can go as high as 6.4 times the Eddington rate.
There are two key pieces of evidence that we have not discussed yet that support this diagnosis.The first, most obviously, is the putative detection of pulsations during a short segment of the XMM-Newton observation 0164560901.If real, this is clear evidence for a NS that is accreting at super-Eddington rates.Our detected signal was at a frequency of 3.32 Hz (corresponding to a period of 0.30 s) with frequency derivative of -1.61 × 10 −7 Hz s −1 and pulse fraction of ∼38 per cent.This was detected with a statistical value of  2 3 ∼ 51, corresponding to a detection significance level of > 99.99 per cent.
However, its detection in only one segment of one out of many observations may mean its pulsations are transient.Similar behaviour has been seen in other PULXs (and, indeed, PULX candidates -see Bachetti et al. 2014;Sathyaprakash et al. 2019;Quintin et al. 2021), although the cause of this behaviour remains a puzzle.One possible clue is that the pulse fraction of the PULX M82 X-2 is proportional to its intensity (Bachetti et al. 2020; see also Lewin et al. 1988).
Interestingly, the data in which we detected the pulsation is the highest flux XMM-Newton observation of the source (cf.Table 3), which might be a factor in the pulse being detected; however it does not explain why this only happens in a short segment of the data, given that the source flux did not change significantly during the observation (see e.g.Fig. 7).One possibility to explain this relates to the scenario proposed by King & Lasota (2020), where they suggest that the pulsations of PULXs are only detectable if the NS spin axis and the accretion disc beaming axes are strongly misaligned; otherwise, the pulse fraction would be small (or tend to zero), and so pulsations would not be detectable.Thus, if the angle between the the NS spin and the accretion disc beaming axes could dynamically change with time, perhaps due to spin axis precession on hour timescales similar to that detected in some NSs (e.g.Heyl & Hernquist 2002;Makishima et al. 2014), a pulse on/off phase might be seen in the data.Alternatively, variations in the obscuration towards the pulsing regions caused by the inner disc structure and/or an outflowing wind along the line of sight might cause the pulsations to become visible for a short period of time (e.g.Kosec et al. 2018;Walton et al. 2018a;Barra et al. 2022).
The other plausible explanation is that the pulsation is not real -indeed, by slicing the multiple datasets into many more datasets, and then searching an extensive region of  −  parameter space, we should perhaps expect to see some false positive signals at relatively high significance.A full analysis of the statistics of false detections is beyond the scope of this paper; however we can at least examine whether, given a candidate signal close to 3.3 Hz, how likely a purely stochastic false signal would be in our data.To do this, we simulated 10,000 fake datasets with the same properties -i.e.mean count rate and exposure time -as the real data segment that the pulsation was detected in, using the hendrics script henfake; the fake data were then searched for a pulsating signal using the same method applied to the real data.We find a probability that random fluctuations would produce a false pulsation of frequency ∼3.3 ± 0.25 Hz with a power of ⩾ 50.6 is 0.001, corresponding to a confidence level of 3.0 that the detected pulsation is real.The other key to the believability of such signals is the detection of another signal at a very similar frequency in another segment/observation (cf.Sathyaprakash et al. 2019).This is hinted at in two other segments, but neither signal is statistically robust.We must therefore conclude that this signal marks NGC 2043 XMM4 as a candidate pulsating object, but this remains to be confirmed.It is, however, another hint that this object may host a NS.
Assuming that the detected pulsation is genuine, compared to other super-Eddington NSs with pulsations -the PULXs -this has the fastest pulsations (marginally faster than, but still comparable with, the ∼ 0.4 s period of NGC 7793 P13) with a relatively high pulse fraction at ∼ 38 per cent.Some theoretical studies have explained the observed properties of PULXs in the context of an accreting magnetar (see e.g.Dall'Osso et al. 2015;Eksi et al. 2015;Mushtukov et al. 2015, but see also King & Lasota 2019), in which it would be expected to spin up as material from the donor is accreted over time.This behaviour has also been seen in the long term (∼ 4 years) monitoring campaign of NGC 7793 P13 (Fürst et al. 2021).In fact, calculations by Meng et al. (2022) predicted that 10 4 years ago the source was a highly magnetised NS with an initial magnetic field strength and pulse period of 10 14 G and ∼ 10 s, respectivelysimilar to those of the other confirmed PULXs; as the source evolves with continuous accretion, the magnetic field decays while its spin becomes faster.This would imply that most detected PULXs are relatively young, compared to NGC 7793 P13 as well as (potentially) NGC 2403 XMM4, which exhibit relatively fast spin periods, but lower magnetic fields (see fig. 2 of Meng et al. 2022).However, NGC 1313 X-2 would be very anomalous in this context given that it is surrounded by a bubble nebula that implies an age in excess of 10 5 years for that PULX (Sathyaprakash et al. 2019).
The other potentially important piece of evidence is the possible detection of absorption features in the spectrum from XMM-Newton observation 0870870501, in particular the feature at 1.68 keV with moderate (2.6) significance (for this discussion, we will disregard the minimally significant feature at 0.99 keV).We first caution that this line energy is close to known instrumental features of XMM-Newton17 , however these features should be properly addressed during the analysis.Indeed, it is clear that we did not see this feature in the other XMM-Newton observations, implying that the feature is intrinsic to the source during that one observation (if it is a real feature).
It is interesting to speculate whether the absorption feature could be a CRSF produced by the strong magnetic field of a NS, similar to those reported by Brightman et al. (2018); Walton et al. (2018b).If it were, it would imply that NGC 2403 XMM4 is a magnetised NS, and would implicitly support the detected pulsation as a genuine feature, since it could be generated by a spinning magnetised NS.If we assume that the line absorption is a result of the cyclotron resonance scattering mechanism, we can determine the magnetic field close to the NS surface.Assuming electron scattering,  12 (G) =  CRSF (keV) × (1 +  grav ) / 11.57, where  12 is the magnetic field in unit of 10 12 G;  CRSF is the line energy in keV; and  grav is gravitational redshift which is assumed to be ∼0.25 (Brightman et al. 2018;Walton et al. 2018b).This predicts a NS surface magnetic field of 1.82 × 10 11 G.However, we note that such a strong magnetic field could also actually beam the emitted X-ray photons, resulting in anisotropic emission (see e.g.Mushtukov & Tsygankov 2022).Assuming a nearly face-on, 2.9M ⊙ NS accreting matter at ∼Eddington rate, a beaming factor of ∼1.3 (or higher for a typical less massive and/or sub-Eddington accreting NS) would be required for the source to reach an average luminosity of 5 × 10 38 erg s −1 .Thus, if the CRSF scenario is correct, we cannot completely rule out that this source is a ∼Eddington or sub-Eddington NS accretor.
However, an alternative interpretation of the 1.68 keV feature (which may also explain the very marginal feature at 0.99 keV) could be that the feature is produced by a combination of atomic features, in both absorption and emission, with the absorption features (at least) resulting from material in the fast ( ∼ 0.2) outflowing wind generated by super-Eddington accretion sources, and so imprinted on its observed X-ray spectrum (see e.g.Pinto et al. 2016).At the spectral resolution of the XMM-Newton EPIC detectors individual features are not often resolved, but broader features resulting from a combination of discrete lines are seen in ULX spectra (Middleton et al. 2014(Middleton et al. , 2015b)).In fact, referring to Fig. 3 of Pinto et al. (2016), it can be seen in the XMM-Newton RGS spectrum of the ULX NGC 1313 X-1 that an absorption line (with > 4 detection significance) has been found at a wavelength of ∼7.5Å (1.65 keV) 18 which can be identified as Ne X with a blueshift of ∼0.2.Interestingly, this absorption line is consistent within errors with the detected feature in NGC 2403 XMM4 at 1.68 ± 0.06 keV.The feature therefore presents possible evidence for an outflowing wind similar in nature to those seen in ULXs; if, similar to ULXs, it is driven from the accretion flow by the extreme radiation pressure of a super-Eddington flow, then this again supports the notion of the accretor being a NS.Indeed, it is notable that both scenarios for the absorption feature outlined here support the presence of a NS accreting at super-Eddington rates.

CONCLUSION
Our interest in NGC 2403 XMM4 was piqued by its inclusion in some older ULX catalogues and studies; however in this work we demonstrate that a classification as a ULX based purely on luminosity is questionable for this source, given that it rarely (if ever) 18 A significant absorption feature with this energy is also detected in the RGS spectrum of the ULX NGC 5408 X-1 (see Fig. 6 of the same paper).
reaches the 10 39 erg s −1 threshold.However, despite not making this threshold, we show that the object displays many behaviours and characteristics that are consistent with super-Eddington accretion: a two-thermal-component spectrum in the 0.3-10 keV regime with a possible steep high energy tail above 10 keV; a soft component whose luminosity and temperature vary together in a manner consistent with an advection-dominated (slim) disc, at least when the spectra are in a similar disc-like state; one observation showing a putative absorption feature, which could be interpreted as a CRSF or absorption in an outflowing wind; and a possible transient pulsation at a similar frequency to PULXs.Although individually each of these results might be questionable -continuum spectra tend to be degenerate to multiple models, the absorption feature is only moderately significant, and the pulsation is only detected in one epoch and so requires a second detection for confirmation it is real -together they all point in one direction, which is that NGC 2403 XMM4 is a NS that is accreting at super-Eddington rates.If so, it is the third such source detected in NGC 2403.

Figure 1 .
Figure1.Twenty year light curve of NGC 2403 XMM4 obtained from Swift (blue boxes), Chandra (magenta diamonds), XMM-Newton (green circles) and NuSTAR (red triangles) observations.The observed luminosity is reported in the 0.3-10 keV energy band, with Chandra and NuSTAR values extrapolated to that band from their best-fitting parameters.

Figure 2 .
Figure 2. Individual XMM-Newton and merged Chandra spectra of NGC 2403 XMM4 with the best-fitting two component model in each case.Top panels: best fitting model (solid-black line; see Section 3.1.2)overlaying the data.The pn, MOS1 and MOS2 data are plotted in blue, red and green, respectively; the dark blue in the bottom-right plot represents the merged Chandra spectrum.Bottom panels: the fitting residuals showing the difference between the model and the data in multiples of , where this is the 1 uncertainty in the value of each data point.

Figure 3 .
Figure 3. Left panel: the best-fitting two component models of individual XMM-Newton and merged Chandra spectra (solid line) along with the contribution from the BB component (dashed line), showing the trend of soft component evolution with the luminosity.Right panel: the relation between the temperature and the intrinsic flux of both individual model components; the data points towards the lower left plotted with a diamond symbol are those of the BB component, while those in the upper right plotted with triangle are the MCD component.The solid black line shows the best-fitting, linear correlation between logarithms of flux and temperature of the BB component of  ∝  1.49 , while the grey dashed and grey dashed-dot lines outline the  ∝  2 and  ∝  4 correlations, respectively.Note that the plotted data are obtained from fitting with the  H tied, excepting for that of XMM-Newton 0870870201 in which the plotted data are adopted from the individual fit (see text for detail.)

Figure 4 .
Figure 4.The absorption features appearing imprinted on the 0870870501 spectrum.The plot definition and colours are as per Fig. 2.

Figure 5 .
Figure 5.The broadband -i.e., XMM-Newton + NuSTAR -spectra of NGC 2403 XMM4 fitted by the three component model (top panel) for observations of epoch 1 (left) and epoch 2 (right).The middle and bottom panels compare the residuals obtained from the best-fitting two and three component models.The orange and purple points present the NuSTAR FPMA and FPMB data, respectively, while the other colour codes and the plot definition are as per Fig. 2.

Figure 6 .
Figure 6.The PSDs of NGC 2403 XMM4 in the full (top), soft (middle) and hard (bottom) energy bands obtained from XMM-Newton observation 0164560901.The horizontal, red-dashed line represents the Poisson noise level.

Figure 7 .
Figure 7.The light curve of NGC 2403 XMM4 with 500 s bin size obtained from XMM-Newton data Obs.ID 0164560901.The red vertical lines represent the boundaries of bad time intervals.The dashed-green line highlights the GTI in which the pulsating signal was detected.

Figure 8 .
Figure 8. Top panel: the  2 3 statistic as a function of frequency, calculated for the segment of XMM-Newton observation 0164560901 in which a pulsation at ∼3.3 Hz is detected.Bottom panel: the folded pulse profile during the period that the pulsating signal is detected; the horizontal dotted line demarcates the mean count rate.

Figure 9 .
Figure 9.The  - parameter space showing the values of  and  determined for the segment of observation 0164560901 in which the source exhibited a pulsation.

Table 2 .
The best fitting results from the merged Swift spectrum, using single component models External absorption column density. Inner disc temperature or photon index for MCD and power-law models, respectively.  2 and number of degrees of freedom for the best fit. Observed X-ray luminosity in the 0.3-10 keV energy band.

Table 3 .
Best fits to the individual XMM-Newton and merged Chandra spectra, as modelled by a BB plus MCD model

Table 4 .
The best fitting results from the broadband spectra modelled by two and three component models