Discovery of a 2.8 s Pulsar in a 2 Day Orbit High-mass X-Ray Binary Powering the Ultraluminous X-Ray Source ULX-7 in M51

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 falling into 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 backgrounds were filtered out.

Table 1.  The XMM-Newton Observations of M51 ULX-7

Data Set	Start Date	pn Eventsa	Period	PF	${T}_{\mathrm{Obs}}$ b	Off-axis Angle	EPIC (${\rm{\Delta }}t$)c
ObsID	(MJD)	(No.)	(s)	(%)	(ks)	(arcmin)	(s)
0112840201	52,654	1241	⋯	<37	20.9 (19.0)	2.3	pn (0.073)
0212480801	53,552	6140	3.2831(2)c	12(2)	49.2 (47.3)	3.4	pn (0.073)
0303420101	53,875	5771	⋯	<17	54.1 (44.0)	3.5	pn (0.073)
0303420201	53,879	6078	⋯	<16	36.8 (34.9)	3.5	pn (0.073)
0677980701	55,719	1206	⋯	<37	13.3 (11.4)	3.6	pn (0.073)
0677980801	55,723	211	2.8014(7)d	82(17)	13.3 (11.4)	3.6	pn (0.073)
0824450901 (A)	58,251	15,082	2.79812(5)	15(2)	78.0 (74.8)	0.0	pn (0.073)
0830191401 (L)	58,263	1037	⋯	<28e	98.0 (94.7)	0.0	pn (0.073), MOS (0.3)
0830191501 (B)	58,282	12,083	2.7977148(2)f	8(2)	63.0 (59.8)	0.0	pn (0.073), MOS (0.3)
0830191601 (C)	58,284	12,559		10(2)	63.0 (59.8)	0.0	pn (0.073), MOS (0.3)

Notes. Uncertainties in the measurements are reported at 1σ confidence level.

^aEvent numbers refer to the total number of photons within the source extraction region for the EPIC pn. ^bDuration of the observation. The effective EPIC pn detector exposure time is given in parentheses. ^cThe number in parentheses indicates the time resolution of the camera(s) used for the timing analysis. ^dNote that the maximum period variation induced by the orbital motion (not corrected for here) is ${\rm{\Delta }}P\sim 1$ ms. ^eEvents from both EPIC pn and MOS have been used in order to infer this upper limit. ^fFor observations 0830191501 and 0830191601, the timing parameters were calculated together; see the text for details.

Download table as:  ASCIITypeset image

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). The times of arrival (ToAs) of the photons were shifted to the barycenter of the solar system with the SAS task BARYCEN (the Chandra position from Kuntz et al. 2016 was used; ${\rm{R}}.{\rm{A}}.\,=\,{13}^{{\rm{h}}}{30}^{{\rm{m}}}01\buildrel{\rm{s}}\over{.} 02$, $\mathrm{decl}.\,=\,47^\circ 13^{\prime} 43\buildrel{\prime\prime}\over{.} 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 the choice of the region to evaluate the background. We finally settled for a circular region of $35^{\prime\prime}$ for the source events and an annular region of inner and outer radii of $50^{\prime\prime}$ and $70^{\prime\prime}$ centered on the source position for the background. Some faint pointlike 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 Section 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 the spectra rebinned to preserve the intrinsic spectral resolution using the SAS tool specgroup.

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

The threshold of 5000 counts is a conservative value based on the relation that links the number of counts to the minimum detectable signal PF in fast Fourier transforms (FFTs), which is given by $\mathrm{PF}={\left\{\left[\tfrac{{P}_{j}}{2M}\right]\tfrac{4}{0.773{N}_{\gamma }}\tfrac{{\left(\pi j/N\right)}^{2}}{{\sin }^{2}(\pi j/N)}\right\}}^{1/2}$, where P_j is the power in the jth Fourier frequency, N_γ and N are the number of counts and bins in the time series, and M is 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 15 M51 sources for which we searched for signals with PASTA, we found that only M51 ULX-7 shows a relatively strong signal at ∼2.8 s with a PF of ∼10% (data set A in Table 1 and the left panel of Figure 1; the 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 $\dot{P}\sim 9.7\times {10}^{-8}\,{\rm{s}}\,{{\rm{s}}}^{-1}$, for which the signal showed a power of about 80 in the corresponding PSD (see upper left and middle panels of Figure 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{int}}$ but rather to the superposition of ${\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\gt 0^\circ$). We then used our Search for Orbital Periods with Acceleration (SOPA) 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 parameter grid, we obtained a first-order value for the orbital period ${P}_{\mathrm{orb}}\sim 2$ days and a projection of the semimajor axis ${a}_{{\rm{X}}}\sin i\sim 25$ lt-s.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We scanned the parameter space describing a pulsar in a circular orbit around the candidate 2.8 s period signal in the data sets of each single XMM-Newton observation. We used only pn data, unless the observing mode of the MOS had a frame time shorter than the spin period, in which case we combined all EPIC data. We found a strong signal only in the UNSEeN observation and two subsequent DPS observations, taken on 2018 May 13, June 13, and June 15, respectively (observation IDs 0824450901, 0830191501, 0830191601, labeled 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 and C). All of the uncertainties reported in this section have a confidence level of 1σ.

Figure 3 shows the orbital period (${P}_{{\rm{o}}{\rm{r}}{\rm{b}}}=1.9969(7)\,{\rm{d}}{\rm{a}}{\rm{y}}{\rm{s}}$) and the projected semimajor axis of the NS orbit, ${a}_{{\rm{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}_{{\rm{a}}{\rm{s}}{\rm{c}}}=58,267.036(6){\rm{M}}{\rm{J}}{\rm{D}}$. Unfortunately, the signal is not strong enough to take 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\lt 0.22$ at a 2σ 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 semimajor axis (${a}_{{\rm{X}}}\sin i\sim 0.5$ lt-s) and 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 }$ NS and the orbital parameters in Table 2, ${M}_{\star }\gt 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,\tfrac{\pi }{2}]$ interval, the probability of having an angle $i\lesssim 26^\circ$ is only 10%. Therefore, we obtain ${M}_{\star }\lesssim 80{M}_{\odot }$ at the 90% confidence level. Similarly, for the average value of the sine function, we find ${M}_{\star }\simeq 27{M}_{\odot }$.

Zoom In Zoom Out Reset image size

Table 2.  Timing Solutions for M51 ULX-7

Parameters	A	B + C
Epoch T0 (MJD)	58,252.36733583	58,283.4440025
Validity (MJD-58,251)	0.92–1.79	31.09–33.78
$P({T}_{0})$ (s)	2.79812(7)	2.79771475(25)
$| \dot{P}({T}_{0})| \times {10}^{-10}$	<26	−2.4(7)
$\nu ({T}_{0})$ (Hz)	0.357383(9)	0.35743458(3)
$| \dot{\nu }({T}_{0})| \times {10}^{-11}$ (Hz s−1)	<35	3.1(8)
${P}_{\mathrm{orb}}$ (days)	1.9969(7)
${a}_{{\rm{X}}}\sin i$ (lt-s)	28.3(4)
${T}_{\mathrm{asc}}$ (MJD)	58,285.0084(12)
e	<0.22
Mass function (${M}_{\odot }$)	6.1(3)
Min. companion mass (${M}_{\odot }$)	$8.3(3)$

Note. Figures in parentheses represent the uncertainties in the least significant digits and are all at a confidence level of 1σ. The upper limit on the eccentricity (e) is at 2σ. The minimum companion mass is computed for an NS of $1.4\,{M}_{\odot }$.

Download table as:  ASCIITypeset image

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 by considering them as nuisance parameters through profile likelihood (Murphy & Vaart 2000). In the first observation, we obtain ${P}^{\left(a\right)}=2.79812(5){\rm{s}}$ and $| {\dot{P}}^{\left(a\right)}| \lt 3\times {10}^{-9}\,{\rm{s}}\ {{\rm{s}}}^{-1}$ with reference epoch ${T}_{0}^{\left(a\right)}=58,252.36734\,{\rm{M}}{\rm{J}}{\rm{D}}$. In the second set of observations, we obtain ${P}^{\left(b\right)}=2.79771475(19)\,{\rm{s}}$ and ${\dot{P}}^{\left(b\right)}=-2.4(6)\times {10}^{-10}\,{\rm{s}}\ {{\rm{s}}}^{-1}$ with reference epoch ${T}_{0}^{\left(b\right)}\,=58,283.44400\,{\rm{M}}{\rm{J}}{\rm{D}}$.

In order to better characterize 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 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σ 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 $\mathrm{PF}\gt 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 Figure 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While the pulse shape remains unchanged among pointings, within the uncertainties, the PF dependence on the energy is only partially similar to that observed in other PULXs (see Figure 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 observations B and C, making M51 ULX-7 the first PULX showing the time-dependent behavior 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 panel of Figure 2), which is in agreement with the idea that the soft thermal component (diskbb), the only component that is changing among observations A, B, and C, is less pulsed (see Section 2.6). In all three observations (A, B, and C), the pulse profile at high energies ($E\gt 1$ keV) precedes that at lower energies by 0.1–0.2 in phase, similar to the case of NGC 5907 ULX1 (Israel et al. 2017a).

Concerning $\dot{P}$, we note that while the second set of observations (B and 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 1 month apart, we obtain a spin-up rate $\dot{P}=-1.5(1)\times {10}^{-10}\,{\rm{s}}\ {{\rm{s}}}^{-1}$. This is marginally compatible with both the larger spin-up rate measured during the second epoch, ${\dot{P}}^{\left(b\right)}$, and 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 (ObsID 0830191401) between pointings A and B (see also Section 4).

We analyzed all 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 ran the three timing techniques described above to find and study any significant signal. We detected the spin modulation in two observations (in one of the two with marginal significance) taken in 2005 (0212480801, see Figure 5) and 2011 (0677980801), both at 5σ (single trial) and with formal $\dot{P}\sim {10}^{-8}$ and ∼ ${10}^{-7}\,{\rm{s}}\ {{\rm{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σ and 3.3σ 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 (55,723 MJD) data set is likely spurious (highly nonsinusoidal 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 ∼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\,{\rm{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 ∼0.4 s longer than that in 2018. Correspondingly, we obtain an average first period derivative of $\dot{P}\simeq -{10}^{-9}\,{\rm{s}}\ {{\rm{s}}}^{-1}$, a value in the range of $\dot{P}$ observed in the other PULXs. As this was obtained from a long baseline (13 yr), 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}}_{\sec }$.

We first simultaneously fit the average spectra of each of the four 2018 observations with a model commonly adopted for PULX spectra in the 0.3–10 keV range, consisting of a soft component (a multicolor blackbody disk, diskbb; Mitsuda et al. 1984) plus a higher-temperature blackbody (bbodyrad), both absorbed by a total covering column (tbabs).²⁴ Abundances were set to that of Wilms et al. (2000) and the 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-component model. Initially, we left all spectral parameters free to vary, and we obtained a statistically acceptable fit (${\chi }_{\nu }^{2}/$dof $=\,1.05/1\,739$). We noticed that ${N}_{{\rm{H}}}$ was consistent within the uncertainties with a constant value across different epochs and that the fit was insensitive to changes in 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 spectral best fit has ${\chi }_{\nu }^{2}/$ dof $=\,1.06/1744$, the best-fit model parameters are reported in Table 3, and the spectra are shown in Figure 6 (upper left panel).

Zoom In Zoom Out Reset image size

Table 3.  Best Fit to the Average 0.3–10 keV Spectrum Using tbabs(diskbb+bbodyrad)

	tbabs	diskbbody			bbodyrad
Obs.	${N}_{{\rm{H}}}$	${{kT}}_{\mathrm{in}}$	Radiusa,b	Flux${}_{\mathrm{dbb}}$ c	kT	Radiusb	Flux${}_{\mathrm{bb}}$ c	${F}_{{\rm{X}}}$ d	${L}_{{\rm{X}}}$ b	${\chi }^{2}$/dof
	(1020 cm−2)	(keV)	(km)		(keV)	(km)			(1039 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/1\,744$
L		${0.19}_{-0.01}^{+0.01}$	${1673}_{-22}^{+26}$	${0.33}_{-0.02}^{+0.02}$	(1.33 fix.)	⋯	$\lt 0.03$ e	${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}$

Notes. The fluxes and luminosity are also reported (in the 0.3–10 keV band). Uncertainties are at the 1σ confidence level.

^aAssuming an inclination angle of 60° and not correcting for a color factor. ^bAssuming d = 8.58 Mpc. ^c0.3–10 keV unabsorbed fluxes in units of 10⁻¹³ erg cm⁻² s⁻¹. ^d0.3–10 keV absorbed fluxes in units of 10⁻¹³ erg cm⁻² s⁻¹. ^e3σ upper limit.

Download table as:  ASCIITypeset image

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 curves 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 (∼0.4–0.6), and rise/decay of the spin pulse profile (∼0.6–0.8 and 0.2–0.4; see Figure 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⁻²), and we left all other parameters free to vary. We found that the temperature and 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 among 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 Figure 6 (upper right panels and bottom panels). The phase-resolved spectra also show 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.

Table 4.  Best Fit to the Phase-resolved 0.3–10 keV Spectra Obtained Using a tbabs(diskbb+bbodyrad) Model

		tbabs	diskbbody		bbodyrad
Obs.	Phase	${N}_{{\rm{H}}}$	${{kT}}_{\mathrm{in}}$	Norm.	kT	Radiusa	${\chi }_{\nu }^{2}$/dof
		(1020 cm−2)	(keV)		(keV)	(km)
A	Min.	5.9b	${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.9b	${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.9b	${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}$

Notes. Uncertainties are at the 1σ confidence level.

^aAssuming d = 8.58 Mpc. ^bFixed to the average spectral best-fit value.

Download table as:  ASCIITypeset image

Matter accreting onto the NS surface releases an accretion luminosity of ${L}_{\mathrm{acc}}(R)={GM}\dot{M}/R\simeq 0.1\dot{M}\,{c}^{2}$ (where M is the NS mass, $\dot{M}$is the mass accretion rate, G is the gravitational constant, R is the NS radius, and c is 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}_{{\rm{m}}}$ is smaller than the corotation radius ${r}_{\mathrm{co}}={\left(\tfrac{{{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}_{{\rm{m}}}$ (propeller phase), and a lower accretion luminosity ${L}_{\mathrm{acc}}({r}_{{\rm{m}}})={GM}\dot{M}/{r}_{{\rm{m}}}$ is released. By adopting the standard expression for the magnetospheric radius,

$$\begin{eqnarray*}\begin{array}{rcl}{r}_{{\rm{m}}} & = & \displaystyle \frac{\xi {\mu }^{4/7}}{{\dot{M}}^{2/7}{\left(2{GM}\right)}^{1/7}}\\ & = & 3.3\times {10}^{7}\,{\xi }_{0.5}\,{B}_{12}^{4/7}\,{L}_{39}^{-2/7}\,{R}_{6}^{10/7}\,{M}_{1.4}^{1/7}\,{\rm{cm}},\end{array}\end{eqnarray*}$$

where ξ is a coefficient that depends on the specific model of disk–magnetosphere interaction (see Campana et al. 2018 for a recent estimate), μ is the magnetic dipole moment, B₁₂ is the dipolar magnetic field at the magnetic poles in units of 10¹² G, ${M}_{1.4}$ is the NS mass in units of 1.4 M_⊙, ${\xi }_{0.5}$ is the ξ coefficient in units of 0.5, and R₆ is the NS radius in units of 10⁶ 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⁻¹ (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 170{P}^{2/3}{M}_{1.4}^{1/3}{R}_{6}^{-1}$ (Corbet 1996; Campana & Stella 2000; Campana et al. 2001, 2002, 2018; Mushtukov et al. 2015; Tsygankov et al. 2016).

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.²⁶ 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}_{{\rm{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⁻¹ 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 ∼ 5 × 10¹⁵ 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 Figure 8. In the same figure, the dotted–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).

Zoom In Zoom Out Reset image size

If M51 ULX-7 isotropically emitted a maximum luminosity ${L}_{\max ,\mathrm{iso}}$ ∼10⁴⁰ erg s⁻¹, an NS magnetic field of about ∼ 10¹⁴ G would be required. The rightmost double-arrowed vertical segment in Figure 8 shows the factor of ∼30 luminosity range observed from the source, with the black circle representing the average value $\langle {L}_{{\rm{X}}}\rangle$ (inferred from the luminosity values reported in Earnshaw et al. 2016 and 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\lt 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}_{\min ,\mathrm{acc}}=b\,{L}_{\min ,\mathrm{iso}}$, and that the maximum accretion luminosity corresponding to the maximum isotropic luminosity ${L}_{\max ,\mathrm{acc}}=b\,{L}_{\max ,\mathrm{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 ∼ 10¹³ G (see Figure 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 $\leqslant$ 10¹³ 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 can 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}\lt -{10}^{-9}$ s s⁻¹, 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}\lt \dot{M}\ {r}_{\mathrm{co}}^{2}P/I$ (where $I\sim {10}^{45}$ g cm² is the NS moment of inertia). This, in turn, gives $\dot{M}\gt 3\times {10}^{18}$ g s⁻¹ for the time-averaged accretion rate and, equivalently, $\langle {L}_{\mathrm{acc}}\rangle \gt 3\times {10}^{38}$ erg s⁻¹, implying that $b\gt 1/12$ and $B\gt 8\times {10}^{11}$ G. Therefore, we conclude that M51 ULX-7 is likely a moderately beamed X-ray pulsar ($1/12\lt b\lt 1/4$) accreting at up to ∼20 times the Eddington rate and with a magnetic field between $\sim 8\times {10}^{11}$ and $\sim {10}^{13}$ G (see the gray shaded area in Figure 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 a high magnetic field by column accretion (see Israel et al. 2017a).

An alternative interpretation of PULX properties, based on models of disk-accreting BHs 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 (King 2009). The accretion rate onto the NS surface is self-regulated close to the Eddington limit, such that magnetic column accretion would work for known PULXs without ever resorting to high, magnetar-like B fields (King et al. 2017). We applied the above interpretation by using the prescription in King & Lasota (2019) and the values of the maximum (isotropic) luminosity, spin, and secular spin-up rate measured for M51 ULX-7. We derived the 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 HMXBs (e.g., the cases of Her X-1, SMC X-1, LMC X-4, etc.; Hickox et al. 2004). Its luminosity, which 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 ∼5–10 times larger than the NS radius, may itself originate 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 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 toward $p\sim 0.5$. Hence, from our spectral analysis, we cannot exclude the existence of a thick disk in this system.

The observed luminosity swing—also encompassing the disappearance (nondetection) of the harder component in the faint state, observation L—together with the suggested ≈40 days of superorbital flux variations (from Swift monitoring; M. Brightman et al. 2020, in preparation), 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 ∼40 day 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, the white region between $b\sim 1/4$ and 1 in Figure 8 would also be allowed, and the NS magnetic field may be as high as ∼ 10¹⁴ 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 ∝${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 ∼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}_{{\rm{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 NS (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.

The timing parameters of M51 ULX-7 firmly place the hosting binary system in the HMXB class, with a companion star of minimum mass ~8 M_☉. 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_☉ and an orbital period of 2 days. In this scenario, the radius and intrinsic luminosity of the donor at this stage are ∼7 ${R}_{\odot }$ and ∼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}}\gt 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 ≈50,000 yr. The estimated average mass transfer rate is largely in excess of the Eddington limit (≳20,000). Assuming a standard accretion efficiency $\eta =({GM}/{c}^{2}R)\sim 10 \%$ for an NS and negligible beaming, the observed X-ray luminosity (∼ 10⁴⁰ erg s⁻¹) implies that only ∼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 disk, 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 ∼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 a relatively large donor-to-accretor mass ratio, q < 15, and with a similar average accreted mass (∼3%; El Mellah et al. 2019; see also Fragos et al. 2015 for an alternative possible scenario).

Source M51 ULX-7 dwells in the same region of 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 RLOF 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 Section 4.1 and for the pulsation dropout in Section 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 days, 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 Section 3 have magnitudes 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 RLOF onto an NS with a donor close to the minimum estimated mass of 8 ${M}_{\odot }$, 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.