Discovery of a 2.8 s pulsar in a 2 d orbit High-Mass X-ray Binary powering the Ultraluminous X-ray source ULX-7 in M51

TL;DR

We discovered 2.8 s pulsations in the ULX-7 source in M51, revealing a high-mass X-ray binary with a neutron star exhibiting secular spin-up and moderate magnetic field, contributing to understanding ultraluminous X-ray sources.

Contribution

First detection of pulsations in ULX-7, establishing it as a high-mass X-ray binary with a neutron star, and providing detailed measurements of its orbital and spin properties.

Findings

ULX-7 hosts a neutron star with a 2.8 s pulsation period.

The system is a high-mass X-ray binary with an OB giant companion.

The neutron star shows secular spin-up consistent with accretion.

Abstract

We discovered 2.8 s pulsations in the X-ray emission of the ultraluminous X-ray source (ULX) M51 ULX-7 within the UNSEeN project, which was designed to hunt for new pulsating ULXs (PULXs) with XMM-Newton. The pulse shape is sinusoidal and large variations of its amplitude were observed even within single exposures (pulsed fraction from less than 5% to 20%). M51 ULX-7 is a variable source, generally observed at an X-ray luminosity between $10^{39}$ and $10^{40}$ erg s$^{-1}$, located in the outskirts of the spiral galaxy M51a at a distance of 8.6 Mpc. According to our analysis, the X-ray pulsar orbits in a 2-d binary with a projected semi-major axis $a_\mathrm{X} \sin i \simeq$ 28 lt-s. For a neutron star (NS) of 1.4 $M_{\odot}$, this implies a lower limit on the companion mass of 8 $M_{\odot}$, placing the system hosting M51 ULX-7 in the high-mass X-ray binary class. The barycentric pulse period decreased by $\simeq$0.4 ms in the 31 d spanned by our May -- June 2018 observations, corresponding to a spin-up rate $\dot{P} \simeq -1.5\times10^{-10}\text{s s}^{-1}$. In an archival 2005 XMM-Newton exposure, we measured a spin period of $\sim$3.3 s, indicating a secular spin-up of $\dot{P}_{\mathrm{sec}}\simeq -10^{-9}\text{ s s}^{-1}$, a value in the range of other known PULXs. Our findings suggest that the system consists of an OB giant and a moderately magnetic (dipole field component in the range $10^{12}$ G $\lesssim B_{\mathrm{dip}}\lesssim 10^{13}$G) accreting NS with weakly beamed emission ($1/12\lesssim b\lesssim1/4$).

Tables

Table 1. Table 1: The XMM–Newton  observations of M51 ULX-7 . Uncertainties in the measurements are reported at 1 σ 𝜎 \sigma confidence level.

Data set	Start date	pn Eventsaafootnotemark:	Period	PF	$T_{\mathrm{Obs}}$bbfootnotemark:	Off-axis angle	EPIC ($\Delta t$)ccfootnotemark:
Obs.ID	(MJD)	(#)	(s)	(%)	(ks)	(arc min)	(s)
0112840201	52654	1241	–	37	20.9 (19.0)	2.3	pn (0.073)
0212480801	53552	6140	3.2831(2)ddfootnotemark:	12(2)	49.2 (47.3)	3.4	pn (0.073)
0303420101	53875	5771	–	17	54.1 (44.0)	3.5	pn (0.073)
0303420201	53879	6078	–	16	36.8 (34.9)	3.5	pn (0.073)
0677980701	55719	1206	–	37	13.3 (11.4)	3.6	pn (0.073)
0677980801	55723	211	2.8014(7)ddfootnotemark:	82(17)	13.3 (11.4)	3.6	pn (0.073)
0824450901 (A)	58251	15082	2.79812(5)	15(2)	78.0 (74.8)	0.0	pn (0.073)
0830191401 (L)	58263	1037	–	28eefootnotemark:	98.0 (94.7)	0.0	pn (0.073), MOS (0.3)
0830191501 (B)	58282	12083	2.7977148(2)fffootnotemark:	8(2)	63.0 (59.8)	0.0	pn (0.073), MOS (0.3)
0830191601 (C)	58284	12559		10(2)	63.0 (59.8)	0.0	pn (0.073), MOS (0.3)

Table 2. Table 2: Timing solutions for M51 ULX-7 .

Parameters	A	B + C
Epoch $T_0$ (MJD)	58252.36733583	58283.4440025
Validity (MJD-58251)	0.92–1.79	31.09–33.78
$P\left(T_0\right)$(s)	2.79812(7)	2.79771475(25)
$\left|\dot{P}\left(T_0\right)\right|\times {10}^{-10}$	26	$-$2.4(7)
$\nu \left(T_0\right)$ (Hz)	0.357383(9)	0.35743458(3)
$\left|\dot{\nu }\left(T_0\right)\right|\times {10}^{-11}$ (Hz s-1)	35	3.1(8)
$P_{\mathrm{orb}}$ (d)	1.9969(7)
$a_X\sin i$ (lt-s)	28.3(4)
$T_{\mathrm{asc}}$ (MJD)	58285.0084(12)
$e$	0.22
Mass function ($M_{\odot }$)	6.1(3)
Min. companion mass ($M_{\odot }$)	$8.3\left(3\right)$

Table 3. Table 3: Best fit to the average 0.3–10 keV spectrum using tbabs(diskbb + bbodyrad) . The fluxes and the luminosity are also reported (in the 0.3–10 keV band). Uncertainties are at 1 ​ σ 1 𝜎 1\sigma confidence level.

	tbabs	diskbbody			bbodyrad
Obs.	$N_H$	$k{T}_{\mathrm{in}}$	Radiusa,b	Flux${}_{}{}^{\ast }{}_{\mathrm{dbb}}{}^{}$	$kT$	Radiusb	Flux${}_{}{}^{\ast }{}_{\mathrm{bb}}{}^{}$	$F_X^{+}$	$L_X$ b	${\chi }^2$/dof
	(${10}^{20}$ cm-2)	(keV)	(km)		(keV)	(km)			(${10}^{39}$ erg s-1)
A	${5.9}_{-0.5}^{+0.6}$	${0.40}_{-0.01}^{+0.01}$	${867}_{-58}^{+58}$	${2.34}_{-0.05}^{+0.05}$	${1.33}_{-0.03}^{+0.04}$	${96}_{-4}^{+4}$	${4.01}_{-0.08}^{+0.08}$	${5.9}_{-0.1}^{+0.1}$	${5.6}_{-0.1}^{+0.1}$	$1.06/\mathrm{1\hspace{0.17em}744}$
L		${0.19}_{-0.01}^{+0.01}$	${1673}_{-22}^{+26}$	${0.33}_{-0.02}^{+0.02}$	(1.33 fix.)	–		${0.21}_{-0.01}^{+0.01}$	${0.29}_{-0.02}^{+0.02}$
B		${0.47}_{-0.02}^{+0.02}$	${728}_{-42}^{+49}$	${3.10}_{-0.07}^{+0.07}$	${1.50}_{-0.05}^{+0.05}$	${86}_{-5}^{+5}$	${4.94}_{-0.12}^{+0.12}$	${7.5}_{-0.1}^{+0.1}$	${7.1}_{-0.1}^{+0.1}$
C		${0.47}_{-0.02}^{+0.02}$	${697}_{-51}^{+51}$	${2.93}_{-0.07}^{+0.07}$	${1.40}_{-0.04}^{+0.05}$	${97}_{-6}^{+6}$	${4.85}_{-0.11}^{+0.11}$	${7.3}_{-0.1}^{+0.1}$	${6.8}_{-0.1}^{+0.1}$

Table 4. Table 4: Best fit to the phase-resolved 0.3–10 keV spectra obtained using a tbabs(diskbb+bbodyrad) model. Uncertainties are at 1 ​ σ 1 𝜎 1\sigma confidence level.

		tbabs	diskbbody		bbodyrad
Obs.	Phase	$N_H$	$k{T}_{\mathrm{in}}$	Norm	$kT$	Radiusa	${\chi }_{\nu }^2$/dof
		(${10}^{20}$ cm-2)	(keV)		(keV)	(km)
A	min.	${5.9}^{\ast }$	${0.41}_{-0.02}^{+0.02}$	${0.45}_{-0.05}^{+0.06}$	${1.43}_{-0.08}^{+0.09}$	${71}_{-6}^{+8}$	$1.07/467$
	raise/decay				${1.27}_{-0.05}^{+0.06}$	${107}_{-8}^{+7}$
	max.				${1.35}_{-0.05}^{+0.05}$	${110}_{-7}^{+7}$
B	min.	${5.9}^{\ast }$	${0.53}_{-0.02}^{+0.02}$	${0.18}_{-0.02}^{+0.02}$	${1.57}_{-0.08}^{+0.09}$	${72}_{-8}^{+8}$	$1.10/437$
	raise/decay				${1.52}_{-0.07}^{+0.07}$	${84}_{-6}^{+6}$
	max.				${1.34}_{-0.06}^{+0.07}$	${104}_{-7}^{+7}$
C	min.	${5.9}^{\ast }$	${0.47}_{-0.01}^{+0.02}$	${0.33}_{-0.04}^{+0.04}$	${1.46}_{-0.07}^{+0.08}$	${81}_{-6}^{+6}$	$1.04/690$
	raise/decay				${1.38}_{-0.05}^{+0.05}$	${99}_{-7}^{+7}$
	max.				${1.29}_{-0.05}^{+0.05}$	${122}_{-8}^{+8}$

Equations

$$r_{\mathrm{m}}$$

Full text

accretion, accretion disks — galaxies: individual: M51a — pulsars: individual (M51 ULX-7, aka CXOM51 J133001.0+471344) — stars: neutron — X-rays: binaries

1 Introduction

Ultraluminous X-ray sources are off-nucleus objects detected in nearby galaxies with X-ray luminosities in excess of $10^{39}$ erg s*-1*, the Eddington luminosity ($L_{\mathrm{Edd}}$) of a 10 $M_{\odot}$ object (e.g. Kaaret et al., 2017, for a review). Under the assumption of a stationary, spherically symmetric accretion flow, $L_{\mathrm{Edd}}$ sets an upper limit to the accretion luminosity that a compact object can steadily produce, since for higher values the accretion flow would be halted by radiation forces. For a compact object accreting fully ionized hydrogen the above limit can be written as $L_{\mathrm{Edd}}=4\pi cGMm_{\mathrm{p}}/\sigma_{\mathrm{T}}\simeq 1.3\times 10^{38}(M/M_{\odot})$ erg s*-1*, where $\sigma_{\mathrm{T}}$ is the Thomson scattering cross section, $m_{\mathrm{p}}$ is the proton mass, and $M$ is the mass of the compact object. Since early discoveries in the ‘70s with the Einstein mission (Long & van Speybroeck, 1983; Fabbiano et al., 1992), the high luminosity of ULXs has been interpreted in terms of accretion at or above the Eddington limit onto BHs of stellar origin ($<$80–100 $M_{\odot}$; e.g. Stobbart et al. 2006; Roberts 2007; Zampieri & Roberts 2009; Feng & Soria 2011), or sub-Eddington accretion onto intermediate-mass BHs (IMBH, $10^{3}$ – $10^{5}$ $M_{\odot}$; e.g. Colbert & Mushotzky 1999; Sutton et al. 2012).

The recent discovery of coherent pulsations with periods between 0.4 and 40 s in the X-ray light curves of four distinct ULXs with luminosities in the $10^{39}$ – $10^{41}$ erg s*-1* range, unambiguously associate these ULXs with accreting NSs, i.e. compact objects with mass of only $\sim$1–2 $M_{\odot}$ (Bachetti et al., 2014; Fürst et al., 2016; Israel et al., 2017a, b; Carpano et al., 2018). Recently a new candidate PULXs, NGC 1313 X-2, has been reported based on weak and transient pulsations at a period of $\sim$1.5 s (Sathyaprakash et al., 2019). These X-ray pulsars demonstrate that accreting NSs can attain extreme luminosities, above 500 times $L_{\mathrm{Edd}}$, which are difficult to interpret in the context of standard accretion models for NS X-ray binaries. A significantly super-Eddington luminosity can be achieved if the magnetic field of the NS is very high, as a result of a marked reduction of the opacities for extraordinary photons: in particular, a luminosity of $\sim$500 $L_{\mathrm{Edd}}$ can be attained for a field strength $>$ $10^{15}$ G (Mushtukov et al., 2015; Dall’Osso et al., 2015), which is expected in magnetars (see e.g. Turolla et al., 2015, for a review). Rotation of the NS and its magnetosphere drags matter at the magnetospheric boundary; if rotation is fast enough the centrifugal force exceeds the gravitational force locally and inhibition of accretion on the NS surface results from the so-called propeller mechanism (Illarionov & Sunyaev, 1975; Stella et al., 1986). Owing to their relatively fast spin period, invoking very strong magnetic fields would imply that the propeller mechanism operates in PULXs; this can be mitigated by assuming that the emission is beamed. In the most luminous PULX, namely NGC 5907 ULX, a beaming factor of 1/100 would be needed (Israel et al., 2017a).

Several scenarios have been proposed to account for the PULX properties. The presence of a strong multipolar magnetic field ($\sim$1014 G) close to the surface of the NSs coupled with a modest degree of beaming appears as a reasonable way out of the problem (Israel et al., 2017a; Chashkina et al., 2017). In this scenario, a standard magnetic dipole field of $\sim$10${}^{12}-$1013 G dominates, at large distances, over the multipolar component, the effect of which are limited to the region close to the surface and to the accretion column base. However, this scenario does not necessarily invoke the presence of magnetars, since magnetars are not merely young NSs with a high dipolar field and/or with a multipolar component; they are also characterized by specific bursting and flaring activity (see e.g. Esposito et al. 2018). Alternatively, standard magnetic fields (in the $10^{11}-10^{13}$ G range, without any significant multipolar component) are envisaged by models in which the disk is fed at a super-Eddington rate, the excess supply is ejected away, and the emission is highly beamed in a geometrically thick inner disk funnel (see e.g. King & Lasota, 2016; King et al., 2017; Pintore et al., 2017; Koliopanos et al., 2017; Walton et al., 2018; king19). We note that two more extragalactic transient pulsators share striking similarities with the above group, i.e. super-Eddington luminosities and/or large first derivative of the spin period, namely XMMU J031747.5–663010 in NGC1313 ($\sim$766 s; Trudolyubov 2008) and CXOU J073709.1+653544 in NGC2403 ($\sim$18 s; Trudolyubov et al. 2007). We suggest here for the first time that they could be two other PULXs that went unnoticed so far.

The discovery of PULXs calls into question the nature ULXs, many of which have been classified as accreting black holes due to their high luminosity: in fact an unknown, but possibly large fraction of ULXs may host an accreting NS rather than a BH (Wiktorowicz et al., 2017; Middleton & King, 2017a, b). Therefore, assessing the nature of the compact objects hosted in ULXs is a key point in understanding the ULX population. In general, the unambiguous identification of the presence of a NS is achieved in X-rays by means of the detection of coherent periodic signals reflecting the spin period of the NS. However, given the small average pulsed fractions (PFs)111Throughout this work, we define the PF as the semi-amplitude of a sinusoidal fit to the pulse profile divided by the source average count rate. of the flux observed in most PULXs (5–15% range), sensitive searches for pulsations require a large number of counts ($\sim$10,000 or higher; see below). For the large majority of the ULXs observed with XMM–Newton, the number of counts collected is by far too low for the detection of pulsations with PFs as small as those observed so far for the fastest-spinning PULXs. Sufficient statistics are currently available only for about $\sim$15 ULXs, out of about 300 known observed with XMM–Newton (Earnshaw et al., 2019); and remarkably, among all the ULXs with good enough statistics data sets $\sim$25% of them are proven to be NSs.

To increase the number of ULXs for which a sensitive search for pulsations can be carried out and to provide a first constraint on the incidence of NSs in ULXs, we observed for long exposures ($\sim 100$ ks) with XMM–Newton the fields of eight nearby galaxies ($3\,\mathrm{Mpc}\leq d\leq 30$ Mpc) hosting a considerable number of ULXs (XMM–Newton Large Project UNSEeN; Ultraluminous Neutron Star Extragalactic populatioN). Among these galaxies is Messier 51a (M51a), also known as NGC 5194 or the Whirlpool Galaxy, a face-on spiral interacting with the dwarf galaxy M51b (NGC 5195). It is located at a distance of $(8.58\pm 0.10)$ Mpc (McQuinn et al., 2016) and hosts a large number of X-ray sources, including nine ULXs (Terashima & Wilson, 2004; Brightman et al., 2018). Based on optical studies, M51a was classified as a Seyfert 2 galaxy (Stauffer, 1982).

M51 ULX-7 (also known as NGC 5194 X-7, Roberts & Warwick 2000; CXOM51 J133001.0+47134, Terashima & Wilson 2004; NGC 5194/5 ULX7 Liu & Mirabel 2005) was first detected with the Einstein Observatory observations at a luminosity above $10^{39}$ erg s*-1* (Palumbo et al., 1985) and located at an offset of about 2.3 arcmin from the central AGN, on the outskirts of a young open cluster on a spiral arm of M51. Deep observations with Chandra showed pronounced variability ($\Delta L_{\mathrm{X}}/L_{\mathrm{X}}\geq 10$) and the presence of a $\sim$7,620 s period modulation (at high $L_{\mathrm{X}}$; Liu et al. 2002). A flux modulation was also observed in an XMM–Newton exposure, though at a significantly different period of $\sim$5,900 s (Dewangan et al., 2005). The variation in period strongly argued against an orbital origin and suggested the presence of some kind of quasi periodic oscillations (QPOs). The XMM–Newton spectral properties and the fact that the source resides near a young massive star cluster with age $T\sim 12$ Myr (Abolmasov et al., 2007) suggested instead that M51 ULX-7 is a high mass X-ray binary (HMXB). More recently, based on a multi-wavelength study (from radio to hard X-rays), it was proposed that the source might be an IMBH accreting in the hard state; however an accreting NS could not be excluded (Earnshaw et al., 2016).

In Section 2 we report on the discovery of coherent pulsations in the X-ray flux of M51 ULX-7 at a period of about 2.8 s with a highly variable amplitude, unambiguously making this source a new member of the rapidly growing class of PULXs. Refined timing analysis allows us to infer an orbital period of about 2 d and a lower limit to the companion star mass of about $8\,M_{\odot}$. Furthermore, the analysis of XMM–Newton archival data makes it possible to infer the secular first period derivative of this new PULX. In Section 3 we make use of the most updated X-ray position of M51 ULX-7 in order to further constrain the optical properties and the nature of its possible counterpart. Finally, in Section 4 we put the inferred properties of M51 ULX-7 in the more general context of the proposed accretion models and of the possible nature of the binary system.

2 XMM–Newton data

2.1 Observations and data reduction

Within the XMM–Newton Large Program UNSEeN, we obtained a 78 ks-long observation of the M51 galaxy followed by three more pointings (one 98 ks and two 63 ks-long) carried out about a month apart as part of the Discretionary time of the Project Scientist (DPS). Before our campaign, the galaxy had already been observed by XMM–Newton on six occasions with shorter exposure times (see Table 1). In the archival observations, the EPIC pn and MOS cameras (Strüder et al., 2001; Turner et al., 2001) were operated in various modes, with different time resolutions (MOS), sizes of the field of view (MOS), and position angles (the UNSEeN field of view was set to avoid targeted sources to fall in CCD gaps). The Science Analysis Software (SAS) v.17.0.0 was used to process the raw observation data files. Intervals of time with anomalously high particle background were filtered out.

For the timing analysis, photon event lists of each source were extracted from circular regions with a radius aimed at minimizing the spurious contribution of nearby objects and diffuse emission in the crowded field of M51. The background was estimated from a nearby, source-free circular region with the same radius. For the timing analysis we mainly used the pn data (with a time resolution of 73 ms), but also data acquired by the MOSs were used when the cameras were operated in Small Window mode (time resolution of 300 ms, that is, only in the three DPS observations). Times of arrival (ToAs) of the photons were shifted to the barycentre of the Solar system with the SAS task BARYCEN (the Chandra position from Kuntz et al. 2016 was used; $\rm RA=13^{h}30^{m}01\fs 02$, $\rm Dec.=47\arcdeg 13\arcmin 43\farcs 8$; J2000).

To study the spectra of M51 ULX-7, both EPIC pn and MOS data were used. Some soft diffuse emission from the host galaxy surrounds the point source. For this reason, we were particularly careful in selecting the source region size and in the choice of the region to evaluate the background. We finally settled for a circular region of 35*′′* for the source events and for an annular region of inner and outer radii of 50*′′* and 70*′′* centered on the source position for the background. Some faint point-like sources lay inside the background region and were excluded from the event selection. We checked that different choices for the region size or the background did not impact the spectral results, in particular in the observation that caught the ULX at the lowest flux (see Sect. 2.6). No significant spectral discrepancies were observed between the pn and MOS spectra inside each observation, therefore we combined the data using the SAS tool epicspeccombine. We checked that the combined spectrum was consistent with the single pn and MOS spectra. The source photons were grouped to a minimum of 25 counts per spectral bin and spectra rebinned to preserve the intrinsic spectral resolution using the SAS tool specgroup.

2.2 Pulsation Discovery

For all time series with at least 5,000 counts, we performed an accelerated search for signals with our PASTA (Pulsation Accelerated Search for Timing Analysis, to be released in a future date) code. PASTA corrects the ToA of each photon accounting for shifts corresponding to period derivatives in the range $-10^{-6}<\dot{P}/P[\mathrm{s}^{-1}]<10^{-6}$, and then looks for peaks above a self-defined detection threshold in the corresponding power spectral density (PSD) even in presence of non-Poissonian noise components (Israel & Stella, 1996).

The threshold of 5,000 counts is a conservative value based on the relation which links the number of counts to the minimum detectable signal PF in fast Fourier transforms (FTTs), which is given by $\mathrm{PF}=\left\{\left[\frac{P_{j}}{2M}\right]\frac{4}{0.773N_{\gamma}}\frac{(\pi j/N)^{2}}{\sin^{2}(\pi j/N)}\right\}^{1/2}$, where $P_{j}$ is the power in the $j$-th Fourier frequency, $N_{\gamma}$ and $N$ the number of counts and bins in the time series, and $M$ the number of averaged FFTs in the final PSD (in our case $M=1$; the formula has been derived from Leahy et al. 1983).

Among the fifteen M51 sources for which we searched for signals with PASTA we found that only M51 ULX-7 shows a relatively strong signal at $\sim$2.8 s with PF of $\sim$10% (dataset A in Table 1 and Fig. 1 left panel; typical PASTA PF upper limits for the other brightest ULXs in M51 during the same observation are in the 7–20% range). PASTA hinted at a $\dot{P}\sim 9.7\times 10^{-8}\text{ s s}^{-1}$ for which the signal showed a power of about 80 in the corresponding PSD (see left and central upper panels of Fig. 2). Based on the fact that all other PULXs are in a binary system, we assume the same is valid for M51 ULX-7. This implies that the observed $\dot{P}$ may not be solely due to the NS intrinsic spin $\dot{P}_{\mathrm{i}nt}$, but instead to the superposition of the $\dot{P}_{\mathrm{int}}$ and an apparent, local $\dot{P}_{\mathrm{orb}}$ caused by the motion of the NS around the center-mass of the binary system (assuming an orbital inclination $i>0\degree$). We then used our SOPA (Search for Orbital Periods with Acceleration) code that performs a similar correction of the ToA as in PASTA but, instead of testing a grid of values for $\dot{P}$, it corrects for a set of values of the orbital parameters, assuming a circular orbit. In a first run with a sparse parameters grid, we obtained a first-order value for the orbital period $P_{\mathrm{orb}}\sim 2$ d and a projection of the semi-major axis $a_{\mathrm{X}}\sin i\sim 25$ lt-s.

2.3 Timing Analysis

We scanned the parameter space describing a pulsar in a circular orbit around the candidate 2.8 s period signal in the datasets of each single XMM–Newton observation. We use only pn data unless the observing mode of the MOS has a frame time shorter than the spin period, in which case we combine all EPIC data. We found a strong signal only in the UNSEeN observation and 2 subsequent DPS observations, taken in May 13, June 13, and June 15, 2018, correspondingly (observation IDs 0824450901, 0830191501, 0830191601, labelled in the following as A, B and C, respectively). Using a direct likelihood technique as described in Israel et al. 2017a (based on Bai 1992 and Cowan et al. 2011), we generated confidence profiles for the orbital parameters in each single observation. Then we combined these results into a single ephemeris, which locks the orbital parameters between observations and allows for two distinct sets of spin parameters ($P$ and $\dot{P}$): one for the first observation in May (A), the other for the pair of observations in June (B, C). All the uncertainties reported in this section have a confidence level of 1$\sigma$.

Figure 3 shows the orbital period ($P_{\mathrm{orb}}=1.9969(7)\text{ d}$) and the projected semi-major axis of the NS orbit, $a_{\mathrm{X}}\sin i=28.3(4)$ lt-s, resulting from this coherent direct likelihood analysis. To complete our description of the circular orbit of the system we estimated the epoch of ascending nodes as $T_{\mathrm{asc}}=58267.036(6)\text{ MJD}$. Unfortunately, the signal is not strong enough to make a further step and probe an elliptical orbit; using the binary model ELL1 (Lange et al., 2001) in TEMPO2 (Hobbs et al., 2006), we set an upper limit on the eccentricity of the orbit of $e<0.22$ at a 2$\sigma$ confidence level. We stress that while the orbital period inferred for M51 ULX-7 is close to that of the revolution of XMM–Newton around the Earth (1.994 days), the projected semi-major axis ($a_{\mathrm{X}}\sin i\sim 0.5$ lt-s) as well as the eccentricity ($e\sim 0.7$) of the spacecraft orbit are so different from those of M51 ULX-7 that we can confidently exclude any relation to it. As we do not see eclipses, the inclination of the system is essentially unconstrained. Figure 3 shows the lower limits on the mass of the companion star for an orbit seen edge-on. Assuming a $1.4\,M_{\odot}$ neutron star and the orbital parameters in Table 2, $\mathbf{M_{\star}>8.3\,M_{\odot}}$. An upper bound on the mass of the companion star can be placed by considering that for a random distribution of orbital inclination angles over the $[0,\frac{\pi}{2}]$ interval, the probability of having an angle $i\lesssim 26\degree$ is only 10%. Therefore, we obtain $M_{\star}\lesssim 80M_{\odot}$ at the 90% confidence level. Similarly, for the average value of the sine function, we find $M_{\star}\simeq\mathbf{27}M_{\odot}$.

If we fix the circular orbit to its best-fit solution, then the spin parameters, $P$ and $\dot{P}$ with a reference epoch at the center of each data set, are strongly constrained, with extremely small uncertainties. However, the spin parameters of both data sets are strongly covariant with the orbital parameters so the stated uncertainties are obtained considering them as nuisance parameters, through profile likelihood (Murphy & Vaart, 2000). In the first observation, we obtain $P^{\left(a\right)}=2.79812(5)\text{ s}$ and $|\dot{P}^{\left(a\right)}|<3\times 10^{-9}\text{ s s}^{-1}$ with reference epoch $T_{0}^{\left(a\right)}=58252.36734\text{ MJD}$. In the second set of observations, we obtain $P^{\left(b\right)}=2.79771475(19)\text{ s}$ and $\dot{P}^{\left(b\right)}=-2.4(6)\times 10^{-10}\text{ s s}^{-1}$ with reference epoch $T_{0}^{\left(b\right)}=58283.44400\text{ MJD}$.

2.4 Signal properties

In order to characterize better the 2.8 s pulsations, we study the time-resolved PFs of the three 2018 observations as a function of the orbital phase.

Figure 4 shows that the PF within each observation is variable and covers a range between about $\sim$18% and 5%. The trend is almost constant during observation A, decreasing in B and increasing in C. Moreover, during the latest part of observation B the signal is not detected with a 3$\sigma$ upper limit of about 3%. This variability seems to be uncorrelated with the orbital phase (signal upper limits in data set B are located at the same orbital phase where pulsations with $\rm PF>10$% are inferred in A and C) suggesting that its origin should not be ascribed to a recurrent (orbital) geometrical effect. The PF distribution for the three data sets suggests that values below 5% are possible, but remain undetected due to low statistics (see the inset in Fig. 4).

While the pulse shape remains unchanged among pointings, within uncertainties, the PFs dependence on the energy is only partially similar to that observed in other PULXs (see Fig. 2 right panels). In fact, the typical increase of the PF for increasing energies is only observed during observation A, while the PFs remain almost constant in observation B and C, making M51 ULX-7 the first PULX showing a time-dependent behaviour of the PF as a function of energy. We also note that the 0.1–12 keV PFs seem to be inversely proportional to the source flux (see Table 1 and right plots of Fig. 2) which is in agreement with the idea that the soft thermal component (diskbb), the only component which is changing among observations A, B and C, is less pulsed (see Sect. 2.6). In all the three observations (A, B and C) the pulse profile at high energies ($E>1$ keV) precedes that at lower energies by 0.1–0.2 in phase, similar to the case of NGC5907 ULX1 (Israel et al., 2017a).

Concerning $\dot{P}$, we note that while the second set of observations (B, C) indicates a spin-up trend, during the first observation (A) the intrinsic $\dot{P}$ could not be constrained. Comparing the spin period at the two epochs, roughly one month apart, we obtain a spin-up rate $\dot{P}=-1.5(1)\times 10^{-10}\text{ s s}^{-1}$. This is marginally compatible with both the larger spin-up rate measured during the second epoch, $\dot{P}^{\left(b\right)}$, and with the $\dot{P}^{\left(a\right)}$ upper limit inferred during the first epoch. This might be a mild indication that the source might have partially entered a propeller spin down phase while it was in a low state during observation L (Obs.ID 0830191401) between pointings A and B (see also Sec. 4).

2.5 Archival data

We analysed all the six XMM–Newton archival data sets of M51 ULX-7 in order to check for the presence of a signal consistent with 2.8 s (see Table 1). In particular, we run the three timing techniques described above to find and to study any significant signal. We detected the spin modulation in two observations (in one of the two with marginal significance) taken in 2005 (021280801) and 2011 (0677980801), both at 5$\sigma$ (single trial) and with formal $\dot{P}\sim 10^{-8}\text{ s s}^{-1}$ and $\sim$ $10^{-7}\text{ s s}^{-1}$, respectively. Taking into account the number of trials, which was estimated as the ratio between the offset in spin frequency from our detections in 2018 and the intrinsic Fourier resolution of each exposure, the significance of the signal goes down to 4.7$\sigma$ and 3.3$\sigma$ in the 2005 and 2011 observations, respectively. An inspection of the timing properties of the candidate signals strongly suggests that the signal detected in the 2011 (55723 MJD) data set is likely spurious (highly non-sinusoidal pulse shape with an abnormally high PF of 82% at variance with the signal properties reported in the above sections; see Table 1). In addition, the 2011 source flux is similar to that observed in deeper observation L, where we did not find any signal down to an upper limit of $\sim$30 % (see Table. 1). Therefore, we considered this detection as not reliable. For the 2005 detection, we cannot disentangle the orbital contribution from the measured spin period $P=3.2831$ s and its time derivative $\dot{P}$, though we can estimate that the maximum $P$ variation induced by the orbital parameters is only of the order of 1 ms, while the 2005 period is $\sim 0.4$ s longer than that in 2018. Correspondingly, we obtain an average first period derivative of $\dot{P}\simeq-10^{-9}\text{ s s}^{-1}$, a value in the range of the $\dot{P}$ observed in the other PULXs. As this was obtained from a long baseline (13 years), it is virtually unaffected by the orbital Doppler shift (which is instead present in each single observation) and, therefore, can be considered as a good estimate of the secular accretion-induced $\dot{P}_{\mathrm{sec}}$.

2.6 Spectral analysis

We first simultaneously fitted the average spectra of each of the four 2018 observations with a model commonly adopted for PULXs spectra in the 0.3-10 keV range, consisting of a soft component (a multicolour blackbody disc; diskbb; Mitsuda et al. 1984) plus a higher-temperature blackbody (bbodyrad), both absorbed by a total covering column (tbabs)222We verified that, in observation B, the source did not spectrally vary during the ‘on’ and ‘off’ phases and therefore we analyzed the whole average source spectrum.. Abundances were set to that of Wilms et al. (2000) and cross section of Verner et al. (1996). We note that the XMM–Newton+NuSTAR spectra of a sample of bright ULXs (comprising two PULXs, Walton et al. 2018) were fitted with the same model plus a cutoffpl component, where the latter was associated with the emission of an accretion column. We have verified that the inclusion of this extra component in our data marginally improved the fit, although its spectral parameters were poorly constrained, because of the smaller energy range. Hence, we proceeded to use only the two-components model. Initially, we let all spectral parameters free to vary and we obtained a statistically acceptable fit ($\chi^{2}_{\nu}/$dof$=1.05/1\,739$). We noticed that the $N_{\rm H}$ was consistent within uncertainties with a constant value across different epochs and that the fit was insensitive to changes of the parameters of the bbodyrad component for the lowest flux observation (L). Therefore, its normalization was set to zero and a new fit was performed by linking $N_{\rm H}$ across all observations (e.g. assuming that the column density did not change significantly between the different epochs). The average spectra best fit has $\chi^{2}_{\nu}/$dof$=1.06/1\,744$, the best fit model parameters are reported in Table 3 and spectra are shown in Fig. 6 (upper left panel).

We subsequently performed a phase-resolved spectral analysis. We selected three phase intervals with respect to the NS spin rotation, where the hardness ratios between the light-curve in the energy bands 0.3–2 and 3–10 keV showed the largest variations (in this case we excluded the time interval of observation B where pulsations are not detected). These phase intervals correspond to the minimum (0.8–1.25), maximum ($\sim$0.4–0.6), and rise/decay of the spin pulse profile ($\sim$0.6–0.8 and 0.2–0.4; see Fig. 2, right panels). For each observation, we fitted the spectra of the three phase-bins simultaneously, using the same model adopted for the average spectra. We fixed the column density to the average spectra best-fit value ($5.9\times 10^{20}$ cm*-2*) and we left all the other parameters free to vary. We found that the temperature and the normalization of the diskbb are constant within each observation with respect to phase changes. The temperature and normalization of the high-energy component are instead variable. For this reason, for each observation, we also linked the diskbb temperature and normalization amongst the spectra of the three spin-phase intervals and performed a further fit. The best-fit values are reported in Table 4 and the corresponding spectra are in Fig. 6 (upper right and bottom panels). The phase-resolved spectra shows also the presence of some broad absorption residuals around 4-5 keV. We tentatively included a gabs model to account for an absorption feature, but the improvement in the best-fit was not significant. More investigations and possibly data would be necessary to confirm or disprove their existence.

3 Archival HST data

We searched for a possible optical counterpart to M51 ULX-7 using archival Hubble Space Telescope (HST) Advanced Camera for Surveys (ACS) data from the Legacy Project targeted at M51 and its companion galaxy NGC 5195 (proposal ID 10452, PI S. Beckwith) and the most accurate source position (Kuntz et al., 2016).

We retrieved calibrated, co-added, geometrically-corrected mosaics from the Hubble Legacy Archive333See https://hla.stsci.edu/hlaview.html. in the F435W, F555W, and F814W filters, corresponding to the B, V, and I bands. We first improved the absolute astrometry of HST images using 445 sources with a match in Gaia Data Release 2 (with an r.m.s. accuracy of $\sim$60 mas). Then, we adopted the X-ray source catalogue produced by Kuntz et al. (2016). To register Chandra astrometry to the HST reference frame, we searched for close, non-ambiguous matches between X-ray sources and HST sources in the F814W band. We selected 19 good reference objects, including 2 bright foreground stars (a further, matching foreground star was excluded from the list because of its large proper motion, measured by Gaia). We then aligned Chandra coordinates so to match that of HST, with a resulting r.m.s. of $\sim 0\farcs 12$, in agreement with results reported by Kuntz et al. (2016). In Figure 7 we show the Chandra error circle of M51 ULX-7 on the HST F555W image, assuming a radius of $0\farcs 38$ (3 times the r.m.s. of the Chandra-to-HST frame registration). Our error circle is broadly consistent with the one adopted by Earnshaw et al. (2016). At least 5 HST sources lie within our error region. The bright object marked as source A, not considered by Earnshaw et al. (2016) as a possible counterpart to M51 ULX-7, has magnitudes $m_{435W}=22.4\pm 0.1$, $m_{555W}=22.4\pm 0.1$, and $m_{814W}=22.5\pm 0.1$, as estimated using the SExtractor software (Bertin & Arnouts, 1996, but a word of caution is appropriate in using photometric results, in view of the very crowded nature of the field). Adopting a distance to M51 of 8.58 Mpc as well as an $E(B-V)=0.030$ (following Earnshaw et al. 2016), its absolute magnitude is $M_{V}\sim-7.4$, with $B-V\sim 0$. Sources 2, 4, 5, and 8 were already discussed by Earnshaw et al. (2016, see their table 6 for magnitudes and colours).

4 Discussion

In 2018, taking advantage of the high throughput of XMM–Newton and good time resolution of the EPIC cameras, within the UNSEeN Large Program, we detected pulsations at a period of $\sim$2.8 s in the X-ray flux of the variable source M51 ULX-7. We found that the signal corresponding to the NS spin is affected by both a secular, intrinsic, spin–up evolution and a Doppler effect due to the revolution of the pulsar around the barycenter of the binary system in an $\sim$2 d-long orbit. All these findings unambiguously point to M51 ULX-7 as a new member of the growing class of PULXs.

During our observations and earlier campaigns (e.g. Earnshaw et al. 2016), the (isotropic) X-ray luminosity of M51 ULX-7 extended over a range from $L_{\mathrm{iso,min}}\leq 3\times 10^{38}$ to $L_{\mathrm{iso,max}}\sim 10^{40}$ erg s*-1*, with variations occurring on timescales longer than several days and an average luminosity level of $\langle L_{\rm iso}\rangle\simeq 4\times 10^{39}$ erg s*-1*.

We detected periodic pulsations in four distinct XMM–Newton observations, three times during May–June 2018 and once in a 2005 archival data set, when the X-ray source luminosity was $L_{\mathrm{X}}\simeq(6$–$8)\times 10^{39}$ erg s*-1*. In two cases (observations B and C), the signal was weaker and could be detected only after photon arrival times were corrected for orbital parameters using our SOPA code described in Sect. 2.2. The pulsed signal showed variable properties that appear to be independent of the orbital phase. In particular, the PF of the spin modulation decreased to 5% within few hours from an initial value of about 12%, and became undetectable (upper limit of about 3%) close to the end of observation B. This is the first time that such strong changes in pulsation amplitude are observed on timescales as short as hours for a PULX. We note that a similar pulsation drop-out has been recently found in NuSTAR data of LMC X-4 when the source was close to the Eddington luminosity (Brumback et al., 2018). During observations B and C, the PF did not show the characteristic increase with energy that was observed in A and in other PULXs, while the B and C X-ray spectra did not show significant variations above a few keV with respect to that of observation A.

The XMM–Newton spectrum of M51 ULX-7 could be modelled with the sum of two thermal components, a commonly adopted model for ULX spectra (see e.g. Pintore et al. 2015; Koliopanos et al. 2017; Walton et al. 2018). We found that the spectral properties during observations A, B and C did not change significantly (spectral parameters broadly consistent to within 3$\sigma$) with the exception of the overall normalization; the unabsorbed 0.3–10 keV luminosity was $\sim$(6–8)$\times$1039 erg s*-1* (for a distance of 8.58 Mpc). The temperature of the soft component, a multicolour blackbody disc, was $\sim$0.4–0.5 keV, in agreement with the soft temperatures observed in other ULX spectra (e.g. Gladstone & Roberts 2009; Sutton et al. 2013; Pintore et al. 2014; Middleton et al. 2015). The soft component may represent the emission from the inner regions of the accretion disc. The large spin-up experienced by the NS implies that the magnetic field truncates the disc at the magnetospheric radius ($r_{\mathrm{m}}\propto B^{4/7}L_{\mathrm{X}}^{-2/7}$, where $B$ is the dipole magnetic field and $L_{\mathrm{X}}$ is the X-ray luminosity). In this scenario, the inner disc radius $R_{\mathrm{in}}$ is equal to $r_{\mathrm{m}}$, implying that $L_{\mathrm{X}}\propto K^{-7/4}$, where $K\propto R_{\mathrm{in}}^{2}$ is the normalization of the diskbb model. We fitted the best-fit $K$ and $L_{\mathrm{X}}$ values with a power law and obtained an index of $-1.6\pm 0.4$, which is fully consistent with the expected value of $-7/4$. The inner disk radius determined from the multicolour blackbody disc ranges between 700 and 1700 km for an assumed disk inclination of 60$\arcdeg$ (higher system inclinations are unlikely owing to the absence of eclipses or dips in the X-ray light-curve), or between 2000 and 5000 applying a colour correction factor with a typical value for this parameter $f=1.7$ (e.g. Miller et al., 2003), that increases the size of the estimated inner disc radius by a factor $f^{2}$. The harder component, which we fitted with a blackbody, had a peak temperature of $\sim$1.4 keV and emitting radius of $\sim$90–100 km; the luminosity associated to it was about 1.6–1.7 times that of the soft component.

The phase-resolved spectroscopy shows that the pulse variability is mainly associated with the harder component, which also drives the energy dependence of the PF. The normalization and temperature of the blackbody change with the pulse phase: the equivalent emitting radius of this component is $\sim$70 and $\sim$110–120 km around the minimum and maximum of the pulse modulation, respectively.

4.1 Accretion model

Matter accreting onto the neutron star surface releases an accretion luminosity of $L_{\rm acc}(R)=GM\dot{M}/R\simeq 0.1\dot{M}\,c^{2}$ (where $M$ is the NS mass, $\dot{M}$ the mass accretion rate, $G$ the gravitational constant, $R$ the NS radius, and $c$ the light speed). According to the standard scenario for accreting, spinning, magnetized NSs, the accreting matter is able to reach the surface of the compact object (and hence produce pulsations at the spin period), when the gravitational force exceeds the centrifugal force caused by drag at the magnetospheric boundary. The latter is estimated based on the magnetic dipole field component, which at large distances from the NS dominates over higher order multipoles. This condition translates into the requirement that the magnetosphere boundary of radius $r_{\mathrm{m}}$ is smaller than the corotation radius $r_{\mathrm{co}}=\left(\frac{GMP^{2}}{4\pi^{2}}\right)^{1/3}$ where a test particle in a Keplerian circular orbit corotates with the central object. When the above condition is not satisfied, accretion is inhibited by the centrifugal barrier at $r_{\mathrm{m}}$ (propeller phase) and a lower accretion luminosity $L_{\mathrm{acc}}{(r_{\mathrm{m}})=GM\dot{M}/r_{\mathrm{m}}}$ is released. By adopting the standard expression for the magnetospheric radius:

[TABLE]

where $\xi$ is a coefficient that depends on the specific model of disk-magnetosphere interaction (see Campana et al. 2018 for a recent estimate), $\mu$ is the magnetic dipole moment , $B_{12}$ is the dipolar magnetic field at the magnetic poles in units of $10^{12}$ G, $M_{1.4}$ the NS mass in units of 1.4 $M_{\odot}$, $\xi_{0.5}$ is the $\xi$ coefficient in units of 0.5 and $R_{6}$ the NS radius in units of $10^{6}$ cm), the minimum accretion luminosity below which the centrifugal barrier begins to operate reads $L_{\mathrm{accr,min}}\simeq 4\times 10^{37}\xi^{7/2}B_{12}^{2}P^{-7/3}M_{1.4}^{-2/3}R_{6}^{5}$ erg s*-1* (here the spin period $P$ is in seconds). The accretion luminosity drop when the NS has fully entered the propeller regime is given by \sim$$170P^{2/3}M_{1.4}^{1/3}R_{6}^{-1} (Corbet, 1996; Campana & Stella, 2000; Campana et al., 2001, 2002; Mushtukov et al., 2015; Tsygankov et al., 2016; Campana et al., 2018).

For the 2.8 s spin period of M51 ULX-7, the accretion luminosity ratio across the accretor/propeller transition is expected to be a factor of about 340, i.e. about 10 times larger than the observed luminosity swing of $L_{\mathrm{iso,max}}/L_{\mathrm{iso,min}}\sim 30$ in the 0.3–10 keV range.444We note that despite careful selection of the source and background regions, see Sect. 2.1, we cannot completely exclude some contamination in our spectra from the diffuse emission of the host galaxy. It is therefore possible at $L_{\mathrm{min}}$ (observation L) the source luminosity was slightly lower than that reported in Table 3. For the luminosity levels at which pulsations are detected, accretion must be taking place uninhibited onto the NS surface. One possibility is that for a fraction of the luminosity range of $L_{\mathrm{iso,max}}/L_{\mathrm{iso,min}}\sim 30$, the source has partially entered the centrifugal gap, with only a fraction of the accretion flow reaching the NS surface and the rest being stopped at $r_{\mathrm{m}}$ (see the case of 4U 0115+63; Campana et al. 2001). In that case, the NS would be close to its equilibrium period, with the magnetospheric radius close to the corotation radius, and undergo large variations in $\dot{P}$ and spin-down episodes (see for example the case of the $\dot{P}$ variation of M82 X-1; Dall’Osso et al. 2015), such that a very high secular spin-up rate would not be expected. This is at variance with the secular $\dot{P}\sim-10^{-9}$ s s*-1* of M51 ULX-7. Consequently, it is likely that the source was in the accretion regime over the entire luminosity swing so far observed.

Calculations by Mushtukov et al. (2015) show that a magnetically-funneled column accretion onto the NS poles can attain highly super-Eddington luminosities for very high magnetic fields (up to a few $10^{3}\,{L_{\mathrm{Edd}}}$ for $B\sim 5\times 10^{15}$ G). The maximum value that can be reached as a function of the NS field at the magnetic poles is plotted as a solid line in Fig. 8. In the same figure, the dot-dashed line separates the region of the accretion regime (on the left) from that of the propeller regime (on the right) for a 2.8 s-spinning accreting NS (our discussion here parallels the one in Israel et al. 2017a).

If M51 ULX-7 emitted isotropically a maximum luminosity $L_{\text{max,iso}}$ $\sim$1040 erg s*-1*, then a NS magnetic field of about $\sim$ $10^{14}$ G would be required. The rightmost double-arrowed vertical segment in Fig. 8 shows the factor of $\sim$30 luminosity range observed from the source, with the black circle representing the averaged value $\langle L_{\rm X}\rangle$ (inferred from the luminosity values reported in Earnshaw et al. 2016 and in this work). It is apparent that for luminosities below $\langle L_{\rm X}\rangle$, the source would straddle the transition to the propeller regime: based on the discussion above, we deem this unlikely.

We consider the possibility that the emission of M51 ULX-7 is collimated over a fraction $b<1$ of the sky. The measured luminosity thus corresponds to the apparent isotropic equivalent luminosity, and the accretion luminosity is thus reduced according to $L_{\mathrm{acc}}=b\ L_{\mathrm{iso}}$. By requiring that the propeller regime has not yet set in at the accretion luminosity corresponding to the minimum detected isotropic luminosity, i.e. $L_{\mathrm{min,acc}}=b\,L_{\mathrm{min,iso}}$, and that maximum accretion luminosity corresponding the maximum isotropic luminosity $L_{\mathrm{max,acc}}=b\,L_{\mathrm{max,iso}}$ is consistent with being produced by column accretion in accordance with Mushtukov et al. (2015), we derive a maximum beaming factor of $b\sim 1/4$ and a maximum NS dipolar magnetic field of $\sim$ $10^{13}$ G (see Fig. 8). For any value of $b$ smaller than $\sim 1/4$ (i.e., a more pronounced degree of beaming) a range of values of $B$ $\leq$1013 G can be found such that the above two requirements are also verified. This may suggest that the beaming factor $b$ can attain very small values and the corresponding accretion luminosity be reduced at will. However an additional constraint comes from the observed spin-up rate which must be sustained by a sufficiently high accretion rate. To work out the minimum accretion luminosity (and hence the accretion rate) required to give rise to a secular spin-up of $\dot{P}<-10^{-9}$ s s*-1*, we consider that the highest specific angular momentum that can be transferred to the NS is that of disk matter entering the NS magnetosphere at the corotation radius $r_{\mathrm{co}}$. The resulting upper limit on the accretion torque translates into the condition $-\dot{P}<\dot{M}\ r_{\mathrm{co}}^{2}P/I$ (where $I\sim 10^{45}$ g cm2 is the NS moment of inertia). This in turn gives $\dot{M}>3\times 10^{18}$ g s*-1* for the time-averaged accretion rate and, equivalently, $\langle L_{\rm acc}\rangle>3\times 10^{38}$ erg s*-1*, implying that $b>1/12$ and $B>8\times 10^{11}$ G. Therefore we conclude that M51 ULX-7 is likely a moderately beamed X-ray pulsar ($1/12<b<1/4$), accreting at up to $\sim$20 times the Eddington rate and with a magnetic field between \sim$$8\times 10^{11} and \sim$$10^{13} G (see the grey-shaded area in Fig. 8). Note that we implicitly assumed that the B-field is purely dipolar. Based on the above interpretation the properties of M51 ULX-7 do not require (but cannot exclude) the presence of higher multipolar components dominating in the vicinity of the star surface; this is unlike the case of the PULX in NGC 5907, whose much larger spin up rate requires a higher $L_{acc}$ and thus powering by column accretion a high magnetic field (see Israel et al. 2017a).

An alternative interpretation of PULX properties, based on models of disk-accreting black holes in ULXs, envisages that a large fraction of the super-Eddington mass inflow through the disk is ejected from within the radius where the disk becomes geometrically thick, with radiation escaping from a collimated funnel (King09). The accretion rate onto the NS surface is self-regulated close to the Eddington limit, such that magnetic-column accretion would work for known PULX without ever resorting to high, magnetar-like B-fields (King et al., 2017). We applied the above interpretation by using the prescription in king19 and the values of the maximum (isotropic) luminosity, spin, secular spin-up rate measured for M51 ULX-7. We derived values of the mass inflow rate through the magnetospheric boundary of about 12 times the Eddington rate and surface magnetic field of $\sim 1.4\times 10^{12}$ G. These values are not in tension with the range of values inferred in our interpretation of M51 ULX-7, despite clear differences in the modeling.

The soft X-ray spectral component of M51 ULX-7 may originate in the accretion disk, as originally proposed to interpret the soft X-ray excess observed in a number of Galactic accreting pulsars in high mass X-ray binaries (e.g. the case of Her X-1, SMC X-1, LMC X-4, etc.; Hickox et al. 2004). Its luminosity, that we estimated from the best-fits of the average spectra of the bright states to be more than 30% of the total, exceeds the energy release of the disk itself by about 2 orders of magnitude; therefore, it is likely powered by reprocessing of the primary central X-ray emission, the hard, pulsed component. Indeed, our diskbb fits to the soft component variations with a color correction (see above)give inner disk radii consistent with the corotation radius of 3300 km. The equivalent blackbody radius of the hard component, being $\sim$5–10 times larger than the NS radius, may originate itself from reprocessing by optically thick curtains of matter in the magnetic funnel that feeds the accretion column from higher altitudes. However, in order to intercept and reprocess $>$30% of the accretion luminosity, the inner disk regions must subtend a large solid angle relative to the central source of primary radiation: it is not clear that such a puffed-up disk would maintain the same surface emissivity law of a standard disk, as implicit in the diskbb model. Substituting the diskbb with a diskpbb model in the fit of the average A, B and C spectra, we could obtain acceptable results with the emissivity index converging towards $p\sim 0.5$. Hence, from our spectral analysis, we cannot exclude the existence of a thick disc in this system.

The observed luminosity swing — encompassing also the disappearance (non detection) of the harder component in the faint state, observation L — together with the suggested $\approx$40 d superorbital flux variations (from a Swift monitoring; Brightman et al., in prep.) may be due to genuine variation of the mass accretion rate onto the NS (as assumed above) or result from partial obscuration of the X-ray emission relative to our line of sight, due, e.g. to a precessing accretion disk (see e.g. Middleton et al. 2018). In the former case a recurrent disk instability or a modulation in the mass transfer rate from the companion induced by a third star in an eccentric $\sim$40 d orbit, might cause the observed luminosity variations. In the latter case the nodal precession of a tilted disk may modulate the observed flux through partial obscuration (see e.g. the so-called slaved, disk model for the superorbital cycle of Her X-1, Roberts 1974). The true luminosity swing of the source would be smaller and confined to the upper range close to $L_{\mathrm{iso,max}}$. In this case also the white region between $b\sim 1/4$ and $b\sim 1$ in Fig. 8 would be allowed and the NS magnetic field may be as high as $\sim$ $10^{14}$ G. In this interpretation, obscuration of the central X-ray source (hard spectral component) may be expected during the lowest flux intervals, as indeed observed in observation L. Being driven by a changing inclination angle of the disk, the apparent luminosity variations of the soft component would be expected to scale approximately with the projected area along the line of sight, such that the inner disk radius from the diskbb fit should be $\propto$ $L^{1/2}$ and its temperature constant. However in our fits the scaling with luminosity is approximately $R_{\mathrm{in}}\propto L^{-0.3}$ and the temperature dropped by a factor of $\sim$2 in observation L. Nevertheless, it might be possible to reproduce these results through a suitably-shaped vertical profile of the precessing disk. Taken at face value the $TR_{\mathrm{in}}\propto L^{-0.3}$ dependence is consistent with the expected dependence of the magnetospheric radius on the (accretion) luminosity, $r_{\mathrm{m}}\propto L^{-2/7}$. This may suggest that the luminosity swing of M51 ULX-7 is driven by changes in the mass accretion rate onto the neutron star (rather than variations of viewing geometry), though in this scenario it is not clear how to account for the disappearance of the hard component in observation L.

4.2 The nature of the system

The timing parameters of M51 ULX-7 firmly place the hosting binary system in the HMXB class, with a companion star of minimum mass $\sim 8\,M_{\odot}$. In the following we assume that the system undergoes Roche lobe overflow (RLOF) with a donor close to the minimum estimated mass of 8$M_{\odot}$ and an orbital period of 2 d. In this scenario, the radius and intrinsic luminosity of the donor at this stage are $\sim$7 $R_{\odot}$ and $\sim$4000 $L_{\odot}$, respectively (e.g. Eggleton 1983; Demircan & Kahraman 1991). However, owing to the large mass ratio ($q=M_{\mathrm{donor}}/M_{\mathrm{NS}}>1$), the mass transfer may be unstable (e.g. Frank et al. 2002), if we approximately consider that the evolution proceeds on a thermal timescale $T_{\mathrm{th}}$, the duration of this phase is $\approx$50,000 yr. The estimated average mass transfer rate is largely in excess of the Eddington limit ($\gtrsim$20,000). Assuming a standard accretion efficiency $\eta=(GM/c^{2}R)\sim 10\%$ for a NS and negligible beaming, the observed X-ray luminosity ($\sim$ $10^{40}$ erg s*-1*) implies that only $\sim$3% of the average transferred mass is in fact accreted onto the NS. Therefore, either the system is close to the propeller stage, and accretion is frequently interrupted/limited, or there are powerful outflows launched from the accretion disc, like those reported recently in a few ULXs (e.g. Pinto et al. 2016; Kosec et al. 2018). For the total duration of the mass transfer phase the total mass deposited onto the NS is $\sim$0.1$M_{\odot}$. These numbers are consistent with the possibility that M51 ULX-7 is an HMXB accreting above the Eddington limit. On the other hand, the evolution of a system like this can become dynamically unstable and rapidly lead to a common envelope phase. Alternatively, a slightly different mechanism, known as wind-RLOF, might be at work in M51 ULX-7, providing a stable mass transfer even for relatively large donor-to-accretor mass ratio, q$<$15, and with a similar average accreted mass ($\sim$3%; elmellah19, see also fragos15 for an alternative possible scenario).

M51 ULX-7 dwells in the same region of the Corbet’s diagram (spin versus orbital periods for accreting pulsars, see Enoto et al. 2014 for a recent compilation) of the OB giant and supergiant HMXB systems in Roche lobe overflow in the Milky Way and the Magellanic Clouds. These objects, Cen X-3 in our Galaxy, SMC X-1 and LMC X-4 (which we mentioned in connection with the superorbital flux and spectral variations in Sec. 4.1 and for the pulsation drop-out in Sec. 2.4), are all rather bright NS pulsators, and the two sources in the Magellanic Clouds, remarkably, shine at, or slightly above, the Eddington luminosity (see e.g. Lutovinov et al. 2013, and Falanga et al. 2015 for their orbital parameters and other characteristics of the systems). We note that the PULX M82 X-2, with a spin period of about 1.4 s and a orbital period of about 2.5 d, lies in the same region of the diagram, though the optical counterpart is unknown. The stellar classification of the donor in M51 ULX-7 is unknown, but some of the candidate counterparts in Earnshaw et al. (2016, figure 12) and in Sect. 3 have magnitude and colors consistent with OB supergiants (typical values are $M_{V}$ from –6 to –7, and $(B-V)$ from –0.1 to –0.4; e.g. Cox 2000).

However, as mentioned above, for Roche lobe overflow onto a NS with a donor close to the minimum estimated mass of 8 M⊙, the system is probably in rapid and unstable evolution, and its optical emission is likely to be far from the expected photometric properties of a single isolated massive star. On the other hand, the optical properties of a wind-RLOF accreting system are not known in detail. In both cases, at the present stage it is difficult to make a meaningful comparison with the observed potential counterpart candidates.

5 Conclusions

M51 ULX-7 is a newly identified PULX in a high mass X-ray binary ($M_{\star}\gtrsim 8M_{\odot}$) with an orbital period of 1.997 days and it is characterized by a spin signal with variable (on time scales of hours) properties. In particular, we note that the detection of such weak (few-percent PF) and variable signals in crowded fields can be challenging, probably limiting our chances to obtain a complete picture of the PULX demography with the current-generation X-ray missions. In this respect, the high throughput and good spatial/timing resolution of Athena are expected to be game changers in determining the NS incidence among ULXs.

Though our observational campaign allowed us to infer a relatively accurate orbital solution, other near future timing observations of M51 ULX-7 are needed to further reduce the uncertainties in the orbital parameters and therefore improve the estimates of the spin parameters in the two epochs considered here. This might firmly establish if the propeller mechanism is at play in this source or not. Besides, a better knowledge of the orbital parameters may turn the marginal detection of pulsation in the two archival observations into robust detections, enabling us to extend back by more than a decade the timing history of this PULX.

Our analysis suggest that a relatively “standard” dipolar magnetic field of 1012–1013 G is sufficient to account for the observed luminosity, though we cannot exclude the presence of a stronger (up to $\sim 10^{14}$ G) multipolar component close to the NS surface.

We thank Dr. Norbert Schartel for granting DPS XMM–Newton time. The scientific results reported in this article are based on observations obtained with XMM–Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. GR acknowledges the support of high performance computing resources awarded by CINECA (MARCONI and GALILEO), under the ISCRA initiative and the INAF-CIENCA MoU; and also the computing centres of INAF – Osservatorio Astronomico di Trieste and Osservatorio Astrofisico di Catania, under the coordination of the CHIPP project, for the availability of computing resources and support. This work has also made use of observations made with the NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive, which is a collaboration between the Space Telescope Science Institute (STScI/NASA), the Space Telescope European Coordinating Facility (ST-ECF/ESA) and the Canadian Astronomy Data Centre (CADC/NRC/CSA). This work has also made use of data from the ESA mission Gaia, processed by the Gaia Data Processing and Analysis Consortium (DPAC). We acknowledge funding in the framework of the project ULTraS ASI–INAF contract N. 2017-14-H.0, project ASI–INAF contract I/037/12/0, and PRIN grant 2017LJ39LM. FB is funded by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska Curie grant agreement no. 664931. FF and CP acknowledge support from ESA Research Fellowships. This research was supported by high performance computing resources at New York University Abu Dhabi. HPE acknowledges support under NASA contract NNG08FD60C. TPR acknowledges support from STFC as part of the consolidated grant ST/K000861/1. We thank the anonymous referee for the helpful comments on the manuscript.