On the compactness of the isolated neutron star RX J0720.4-3125

TL;DR

This study reanalyzed extensive X-ray data of neutron star RX J0720.4-3125 using Bayesian methods, revealing complex, energy-dependent light curves, a revised spin period, and estimates of gravitational redshift and compactness.

Contribution

It introduces a Bayesian approach to analyze neutron star light curves and determines a revised spin period, gravitational redshift, and compactness of RX J0720.4-3125.

Findings

The spin period is likely 16.78 seconds, twice the previously reported value.

Light curves are complex and energy-dependent, best modeled with a double-peaked profile.

Estimated gravitational redshift is approximately 0.205, indicating high compactness.

Abstract

The data from all observations of RX J0720.4-3125 conducted by XMM-Newton EPIC-pn with the same instrumental setup in 2000-2012 were reprocessed to form a homogenous data set of solar barycenter corrected photon arrival times registered from RX J0720.4-3125. A Bayesian method for the search, detection, and estimation of the parameters of an unknown-shaped periodic signal was employed as developed by Gregory & Loredo (1992). A number of complex models (single and double peaked) of light curves from pulsating neutron stars were statistically analyzed. The distribution of phases for the registered photons was calculated by folding the arrival times with the derived spin-period and the resulting distribution of phases approximated with a mixed von Mises distribution, and its parameters were estimated by using the Expected Maximization method. Spin phase-resolved spectra were extracted, and a number of highly magnetized atmosphere models of an INS were used to fit simultaneously, the results were verified via an MCMC approach. The phase-folded light curves in different energy bands with high S/N ratio show a high complexity and variations depending on time and energy. They can be parameterized with a mixed von Mises distribution, i.e. with double-peaked light curve profile showing a dependence of the estimated parameters (mean directions, concentrations, and proportion) upon the energy band, indicating that radiation emerges from at least two emitting areas. The genuine spin-period of the isolated neutron star RX J0720-3125 derived as more likely is twice of that reported in the literature (16.78s instead of 8.39s). The gravitational redshift of RX J0720.4-3125 was determined to $z=0.205_{-0.003}^{+0.006}$ and the compactness was estimated to $(M/M_{Sun})/R(km)=0.105 \pm 0.002 $.

Tables

Table 1. Table 1: XMM- Newton  EPIC-pn Full Frame mode observations of RX J0720.4 − - 3125

Obs. ID	Date begin	MJD	Exposure	Effective	Number	Group
		[start]	[ksec]	exposure [ksec]	of counts	No.
0124100101	2000-05-13 01:42:22	51677.104103	65.9	33.6	205200	I
0156960201	2002-11-06 17:51:49	52584.761786	30.2	19.7	125305	”
0156960401	2002-11-08 19:25:02	52586.826653	32.0	30.0	189113	”
0164560501	2004-05-22 10:15:22	53147.442236	52.0	25.4	166528	II
0300520201	2005-04-28 08:41:05	53488.377978	53.3	30.0	193632	”
0300520301	2005-09-22 23:44:23	53636.004188	53.0	35.2	230808	”
0400140301	2006-05-22 04:44:47	53877.228785	21.9	16.2	103291	III
0400140401	2006-11-05 11:19:29	54044.489372	21.9	20.0	127835	”
0502710201	2007-05-05 17:01:25	54225.732373	21.9	15.0	94214	”
0502710301	2007-11-17 05:14:32	54421.236398	24.9	23.0	142855	”
0554510101	2009-03-21 14:14:24	54911.669656	21.9	14.8	92506	IV
0601170301	2009-09-22 04:27:35	55096.200867	30.8	13.2	81991	”
0650920101	2011-04-10 23:59:26	55662.016898	21.9	16.1	99286	”
0670700201	2011-05-02 23:25:17	55684.067143	28.8	10.1	62887	V
0670700301	2011-10-01 03:47:26	55835.184239	26.9	24.1	150050	”
0690070201	2012-09-18 08:50:53	56188.369444	25.1	24.9	159712	”

Table 2. Table 2: Results of the application of the GL method in the frequency range of 0.05-0.07 Hz for the broad energy band of 0.16-1.60 keV.

Obsid	MJD ∗*∗*$\ast$The mid-time of the pointed observation.	Frequency	HPDlo(68%)	HPDhi(68%)	Most probable
	[mid-time]	(Hz)	[Hz]	[Hz]	number of phase bins
0124100101	51677.485469	0.0595869	0.0595863	0.0595873	12
0156960201	52584.936555	0.0595867	0.0595858	0.0595873	10
0156960401	52587.011838	0.0595868	0.0595859	0.0595875	11
0164560501	53147.743162	0.0595923	0.0595919	0.0595928	12
0300520201	53488.686427	0.0595867	0.0595863	0.0595873	12
0300520301	53636.310901	0.0595867	0.0595860	0.0595872	11
0400140301	53877.355521	0.0595876	0.0595870	0.0595882	12
0400140401	54044.616108	0.0595868	0.0595864	0.0595876	12
0502710201	54225.859109	0.0595874	0.0595851	0.0595878	12
0502710301	54421.380495	0.0595865	0.0595859	0.0595868	10
0554510101	54911.796392	0.0595857	0.0595844	0.0595868	11
0601170301	55096.379108	0.0595869	0.0595858	0.0595876	11
0650920101	55662.143634	0.0595851	0.0595841	0.0595858	4
0670700201	55684.233810	0.0595872	0.0595864	0.0595877	11
0670700301	55835.339910	0.0595873	0.0595866	0.0595877	12
0690070201	56188.514699	0.0595870	0.0595865	0.0595874	10

Table 3. Table 3: Simultaneous fitting results of combined, spin-phase averaged X-ray spectra of RX J0720.4 − - 3125 with spectral model of an absorbed blackbody with and without Gaussian absorption line (in XSPEC tbabs*bbodyrad*gabs )

Fitted	Groups∗*∗*$\ast$See subsection 3.2, Table 1 and Fig 1 for definition of the data groups.
Parameter∗⁣∗**∗⁣∗**$\ast \ast$Parameter uncertainty corresponds to the 90% of the confidence interval. $D_{\mathrm{line}}$ presents the strength of a Gaussian line (for details see XSPEC manual).	I	II	III	IV	V	Remark
Energy range	0.16$-$1.2 keV
$nH\times {10}^{20}\left[{\mathrm{cm}}^{-2}\right]$	$1.06\pm 0.06$	$1.09\pm 0.11$	$0.87\pm 0.10$	$1.02\pm 0.10$	$0.92\pm 0.07$	free
$kT\left[\mathrm{eV}\right]$	$85.2\pm 0.4$	$91.3\pm 0.4$	$91.5\pm 0.6$	$89.5\pm 0.5$	$88.2\pm 0.4$	free
$E_{\mathrm{line}}\left[\mathrm{eV}\right]$	$329.0\pm 3.3$	ibid	ibid	ibid	ibid	linked
${\sigma }_{\mathrm{line}}\left[\mathrm{ev}\right]$	$18.3\pm 8.0$	$70.0\pm 6.7$	$77.2\pm 8.7$	$70.4\pm 8.6$	$69.2\pm 8.4$	free
$D_{\mathrm{line}}$∗⁣∗**∗⁣∗**$\ast \ast$Parameter uncertainty corresponds to the 90% of the confidence interval. $D_{\mathrm{line}}$ presents the strength of a Gaussian line (for details see XSPEC manual).	$0.011\pm 0.002$	$0.076\pm 0.009$	$0.074\pm 0.010$	$0.059\pm 0.096$	$0.044\pm 0.072$	free
Reduced ${\chi }^2$	${\chi }^2/\mathrm{dof}=853.07/581=1.47$
$nH\times {10}^{20}\left[{\mathrm{cm}}^{-2}\right]$	$1.04\pm 0.04$	ibid	ibid	ibid	ibid	linked
$kT\left[\mathrm{eV}\right]$	$85.2\pm 0.4$	$90.6\pm 0.3$	$91.1\pm 0.3$	$89.1\pm 0.3$	$88.0\pm 0.3$	free
$E_{\mathrm{line}}\left[\mathrm{eV}\right]$	$328.0\pm 6.2$	$323.0\pm 3.1$	$340.7\pm 3.7$	$329.9\pm 4.0$	$342.9\pm 4.3$	free
${\sigma }_{\mathrm{line}}\left[\mathrm{ev}\right]$	$26.3\pm 10.9$	$82.7\pm 3.5$	$74.4\pm 5.3$	$76.8\pm 5.3$	$63.7\pm 5.6$	free
$D_{\mathrm{line}}$∗⁣∗**∗⁣∗**$\ast \ast$Parameter uncertainty corresponds to the 90% of the confidence interval. $D_{\mathrm{line}}$ presents the strength of a Gaussian line (for details see XSPEC manual).	$0.012\pm 0.002$	$0.095\pm 0.006$	$0.069\pm 0.007$	$0.066\pm 0.006$	$0.039\pm 0.004$	free
Reduced ${\chi }^2$	${\chi }^2/\mathrm{dof}=835.15/587=1.42$

Table 4. Table 4: Simultaneous fitting results of combined, spin-phase resolved X-ray spectra of RX J0720.4 − - 3125 using a condensed iron surface and hydrogen atmosphere model.

Fitted	Groups∗*∗*$\ast$See subsection 3.2, Table 1, Fig 1, and Fig. 9 for definition of the data groups and parameters.
Parameter∗⁣∗**∗⁣∗**$\ast \ast$Parameter uncertainty corresponds to the 95% of the HPD interval.	I	II	III	IV	V	Remark
$nH\times {10}^{20}\left[{\mathrm{cm}}^{-2}\right]$	$1.20\pm 0.06$	$1.17\pm 0.06$	$1.17\pm 0.06$	$1.18\pm 0.06$	$1.18\pm 0.06$	free
$k{T}_{p1}\left[\mathrm{eV}\right]$	$104.0\pm 4.3$	$110.2\pm 4.1$	$108.3\pm 4.5$	$107.5\pm 4.0$	$106.6\pm 3.6$	free
$k{T}_{p2}\left[\mathrm{eV}\right]$	$102.5\pm 4.2$	$107.7\pm 3.9$	$109.5\pm 4.8$	$105.8\pm 3.8$	$105.6\pm 3.4$	free
$B_{p12}\times {10}^{13}\left[G\right]$	$3.00\pm 0.15$	ibid	ibid	ibid	ibid	linked
$\theta {[}^{°}]$	$88.9\pm 4.8$	ibid	ibid	ibid	ibid	linked
$i{[}^{°}]$	$42.8\pm 2.2$	ibid	ibid	ibid	ibid	linked
$z$	$0.205\pm 0.001$	ibid	ibid	ibid	ibid	linked
${\sigma }_{\mathrm{line}}\left[\mathrm{ev}\right]$	$268.3\pm 50.9$	$313.4\pm 38.6$	$329.0\pm 56.9$	$302.7\pm 60.3$	$310.6\pm 60.4$	free
${\tau }_{\mathrm{line}}$	$0.80\pm 0.24$	$1.56\pm 0.26$	$1.38\pm 0.26$	$1.37\pm 0.26$	$1.18\pm 0.22$	free
Reduced ${\chi }^2$	${\chi }^2/\mathrm{dof}=6070/4910=1.24$
C$-$statistic	6050 using 5050 PHA bins and 4910 dof

Table 5. Table 5: Results of the fitting of spin-phases with a two-component mixed von Mises distribution.

Groups	Energy range	N	${\mu }_1$/${\mu }_2$	${\kappa }_1$/${\kappa }_2$	p
	[kev]	[counts]
GRI	0.16$-$0.5	294516	0.507$\pm$0.003	3.124$\pm$0.367	0.349
”	0.50$-$0.6	30943	0.508$\pm$0.007	2.236$\pm$0.880	0.415
”	0.60$-$0.7	16125	0.479$\pm$0.011	2.032$\pm$0.843	0.408
”	0.70$-$1.2	14338	0.495$\pm$0.013	3.286$\pm$1.714	0.311
GRII	0.16$-$0.5	359528	0.494$\pm$0.002	3.101$\pm$0.282	0.352
”	0.50$-$0.6	51604	0.496$\pm$0.005	2.622$\pm$0.732	0.392
”	0.60$-$0.7	30956	0.514$\pm$0.008	1.333$\pm$0.108	0.511
”	0.70$-$1.2	29568	0.502$\pm$0.006	2.159$\pm$0.483	0.402
GRIII	0.16$-$0.5	234977	0.502$\pm$0.003	3.138$\pm$0.271	0.347
”	0.50$-$0.6	32629	0.500$\pm$0.011	2.397$\pm$1.151	0.381
”	0.60$-$0.7	19045	0.520$\pm$0.009	1.363$\pm$0.294	0.512
”	0.70$-$1.2	17789	0.524$\pm$0.012	1.454$\pm$0.200	0.514
GRIV	0.16$-$0.5	210538	0.508$\pm$0.003	2.796$\pm$0.330	0.360
”	0.50$-$0.6	27712	0.498$\pm$0.008	1.911$\pm$0.728	0.505
”	0.60$-$0.7	15814	0.508$\pm$0.009	1.467$\pm$0.476	0.507
”	0.70$-$1.2	14840	0.508$\pm$0.007	1.346$\pm$0.368	0.506
GRV	0.16$-$0.5	352642	0.497$\pm$0.002	3.130$\pm$0.260	0.342
”	0.50$-$0.6	43187	0.505$\pm$0.005	1.611$\pm$0.461	0.507
”	0.60$-$0.7	23725	0.498$\pm$0.006	1.709$\pm$0.451	0.484
”	0.70$-$1.2	21436	0.508$\pm$0.011	1.439$\pm$0.636	0.507
Student’s t$-$test			H0: ${\mu }_1/{\mu }_2\equiv 0.5$	H0: ${\kappa }_1/{\kappa }_2\equiv 1.0$	H0: $p\equiv 0.5$
P$-$value			0.14	$9.8\times {10}^{-9}$	0.005
Wilcoxon test
P$-$value			0.08	$1.9\times {10}^{-6}$	0.008

Equations

$$p(\nu |D,M_m)=\frac{C}{\nu }{\int }_0^{2\pi }d\varphi \frac{1}{W_m(\nu ,\varphi )},$$

$$A=\frac{{\sum }_i|{\mathrm{CR}}_i-{\mathrm{CR}}_{\mathrm{mean}}|}{{\sum }_i{\mathrm{CR}}_i},$$

Full text

*Methods. *The data from all observations of RX J0720.4$-$3125 conducted by XMM-Newton EPIC-pn with the same instrumental setup in 2000-2012 were reprocessed to form a homogenous data set of solar barycenter corrected photon arrival times registered from RX J0720.4$-$3125. A Bayesian method for the search, detection, and estimation of the parameters of an unknown-shaped periodic signal was employed as developed by Gregory & Loredo (1992). A number of complex models (single and double peaked) of light curves from pulsating neutron stars were statistically analyzed. The distribution of phases for the registered photons was calculated by folding the arrival times with the derived spin-period and the resulting distribution of phases – approximated with a mixed von Mises distribution –, and its parameters were estimated by using the Expected Maximization method. Spin phase-resolved spectra were extracted, and a number of highly magnetized atmosphere models of an INS were used to fit simultaneously, the results were verified via an MCMC approach.

*Results. *The phase-folded light curves in different energy bands with high S/N ratio show a high complexity and variations depending on time and energy. They can be parameterized with a mixed von Mises distribution, i.e. with double-peaked light curve profile showing a dependence of the estimated parameters (mean directions, concentrations, and proportion) upon the energy band, indicating that radiation emerges from at least two emitting areas.

*Conclusions. *The genuine spin-period of the isolated neutron star RX J0720$-$3125 derived as more likely is twice of that reported in the literature (16.78s instead of 8.39s). The gravitational redshift of RX J0720.4$-$3125 was determined to $z=0.205_{-0.003}^{+0.006}$ and the compactness was estimated to $(M/M_{\sun})/(R/\mathrm{km})=0.105\pm 0.002$.

1 Introduction

A comprehensive understanding of neutron stars is one of the key challenges facing astrophysics and astronomy. In particular, the observational determination of the gravitational redshift of an isolated neutron star (INS) can be used to put constraints on the theoretical models of superdense matter. This goal can be achieved, for instance, by modeling and analyzing the observed emergent spectrum of a thermally emitting, pulsating INS via rotational phase-resolved X-ray spectroscopy (Zavlin et al. 1995; Motch et al. 2003; Pérez-Azorín et al. 2006; Zane & Turolla 2006; Ho et al. 2007; Suleimanov et al. 2010; Hambaryan et al. 2011, 2014, also developed and successfully applied by us).

This spectroscopy involves simultaneous fitting of high quality spectra111High S/N ratio (relative errors in each spectral bin are of the order of a few percent) spectra for different rotational phases can be achieved by combining observations performed at different epochs. from different spin phases with models of highly magnetized atmospheres of INSs. This allows us to constrain a number of physical properties of the X-ray emitting areas, including their temperatures, magnetic field strengths at the poles, and their distribution parameters. In addition, we may place some constraints on the viewing geometry of the emerging X-ray emission and the gravitational redshift and, hence, the compactness of the INS. This model of spectra is based on various local models, such as blackbody or condensed (iron or other heavy element) surface, covered by a thin hydrogen layer.

On the other hand, it is clear that for such kind of study the true spin parameter of the rotating INS must be known at different epochs, i.e. X-ray timing analysis is also needed for the application of advanced methods (see Sect. 3.1).

At the beginning of 1990ies, the ROSAT X-ray observatory has discovered a small group (so far 7) nearby, thermally emitting and radio-quiet INSs (e.g. see reviews by Haberl 2007; Mereghetti 2008; Turolla 2009, and references therein). Meanwhile, multi-epoch X-ray observations using the XMM-Newton  and Chandra telescopes allowed to detect X-ray pulsations and to uncover absorption features in the spectra of some of them (Haberl et al. 2003; Hambaryan et al. 2009; Borghese et al. 2015).

RX J0720.4$-$3125 is a special case among the nearby thermally emitting INSs showing long-term variations in its timing and spectral parameters. A clear variation in the spectrum has been detected by the XMM-Newton  high spectral resolution Reflection Grating Spectrometers (RGS; de Vries et al. 2004). This discovery was subsequently confirmed and analyzed by imaging-spectroscopic EPIC-pn observations (Haberl et al. 2006; Hohle et al. 2009).

The spin-phase-averaged X-ray spectra of RX J0720.4$-$3125 can formally be modeled by a blackbody (kT $\sim$84$-$94 eV) plus a broad ($\sim$ 10$-$70 eV wide) absorption feature centered at $\sim$ 300 eV (Haberl et al. 2004; Hohle et al. 2009, 2012b) undergoing rather sensible changes over a timescale of a few years (however, see also Viganò et al. 2014, stating that the observed broad absorption feature, partially, might be owing to the X-ray spectral distortions induced by inhomogeneous temperature distributions of the neutron star surface).

A narrow ($\sim 5$ eV wide) absorption feature was identified in XMM-Newton RGS data at $\sim$ 570 eV, possibly due to highly ionized oxygen of circumstellar origin (Mendez et al. 2004; Hambaryan et al. 2009; Hohle et al. 2012a).

In this paper, we performed a rotational phase-resolved spectral re-analysis of the XMM-Newton EPIC-pn data from the avaliable multi-epoch X-ray observations of RX J0720.4$-$3125.

2 Observations and data reduction

RX J0720.4$-$3125 was observed many times by XMM-Newton with different instrumental setup (Table 1). Here we focus mainly on the high quality data collected with EPIC-pn in Full Frame mode (Strüder et al. 2001) from the 16 XMM-Newton observations with same instrumental setup (positioned on-axis), distributed over 12.4 years, in total presenting about 350 ks of effective exposure time, and consisting of $\sim$ 2.2 Million registered source photons.

The data were reduced using standard threads from the XMM-Newton  data analysis package SAS version 14.0.0. We reprocessed all observations listed in Table 1 with the standard metatask epchain. To determine good time intervals, which were free of background flares, we applied a filter to the background light curves and performed a visual inspection. This reduced the total exposure time by $\sim$ 35%. Solar barycenter corrected source and background photon events files were produced from the cleaned single-pixel events, using an extraction radius of 45$-$60$\arcsec$ depending on the brightness of RX J0720.4$-$3125 in the corresponding pointed observation.

We extracted light curves of RX J0720.4$-$3125  and the corresponding background from nearby, source-free regions. We then used the SAS task epiclccorr to correct the observed count rates for various sorts of detector inefficiencies (vignetting, bad pixels, dead time, effective areas, etc.) in different energy bands for each pointed observation for uninterrupted good time interval exceeding 900 sec (see Fig. 1).

3 Data analysis

In order to extract spin-phase-resolved spectra, one has to determine the rotational phase of each registered photon, i.e. an accurate estimate of the most probable value of the spin-period is required for each epoch.

In the case of RX J0720.4$-$3125, apart from the clear variation detected in the spectrum, a significant variation of the spin-period also was reported. A number of published papers are devoted to the timing analysis and behavior of spectral variations of this enigmatic object (Cropper et al. 2001, 2004; van Kerkwijk et al. 2007; Hohle et al. 2009). To explain significant phase residuals from a steady spin-down model, different interpretations have been suggested; free precession and glitch possibilities are mostly debated (see, e.g. Hohle et al. 2009, 2012a, possible glitch ($\Delta\nu\sim 4.1(12)$ nHz at MJD=52866$\pm$73d)).

However, sometimes an unambiguous identification of the true spin-period of a spinning compact star (where the pulsed emission is dominated from emitting areas with slightly different physical characteristics in comparison to the other parts of the surface or atmosphere, e.g. around the magnetic poles) is very difficult from simple analysis of a periodogram (Buccheri et al. 1983; Ögelman & van den Heuvel 1989)222A commonly used statistics in X-ray astronomy:

$Z_{m}^{2}=\frac{2}{N}{\displaystyle\sum_{k=1}^{m}\left\{\left(\sum_{i=1}^{N}\sin\left(2k\pi\phi_{i}\right)\right)^{2}+\left(\sum_{i=1}^{N}\cos\left(2k\pi\phi_{i}\right)\right)^{2}\right\}}$ with phases $\phi_{i}$ at trial frequency. Here, $\phi_{i}=\nu\,(t_{i}-t_{0})+\dot{\nu}\,(t_{i}-t_{0})^{2}/2$ is the phase, $\nu$ and $\dot{\nu}$ are the trial frequency and its derivative, $t_{i}$ ($i=1,...,N$) are the photon arrival times, and $t_{0}$ is the epoch of zero phase. Note, $Z_{m}^{2}$ is a sufficient statistic for detection of sinusoidal signals in comparison to the arrival times, they are distributed uniformly (obeying constant-rate Poissonian counting statistic, i.e. a periodic component is absent). , i.e. by assigning the frequency corresponding to the most significant peak in a periodogram to the true rotational period.

The matter is that, the genuine spin frequency in a periodogram highly depends on the unknown true light curve shape (e.g. single-, double-, or multiple-peaked profile), which in turn depends on the physical characteristics of the emitting areas (e.g. temperatures, sizes, locations, i.e. ratios and locations of the maxima and minima in a light curve) as well from viewing geometry (angles between rotation, magnetic axes, line of sight, gravitational redshift, etc., e.g., see Poutanen & Beloborodov 2006).

On the other hand, it is clear that if there is a periodic signal in the data set with frequency $\nu$, then there is also a less significant signal at $\mathrm{0.5}\nu$. Therefore, the most significant peak in the periodogram unequivocally corresponds to the genuine spin frequency for a single-peaked light curve shape. The situation is not simple for the case of a signal from a double- or multiple-peaked light curve shape. In this case, depending on the parameters mentioned above, the periodogram can be identical to a single-peaked one. Moreover, the dependence of an applied statistic (e.g. $\mathrm{Z^{2}_{m}}$, see also Appendix A) from a trial frequency will show a highly significant detection only at twice of the true spin frequency ($\mathrm{2}\nu$), if the two emitting areas have identical physical characteristics and favorable viewing geometry (e.g. orthogonal rotator, see expression for $\mathrm{Z^{2}_{1}}$ and Fig. 2 in Hambaryan et al. 2015).

These arguments have motivated us to perform a complete X-ray timing re-analysis of the observational data sets of RX J0720.4$-$3125 for the identification of the true pulsation period by using a more advanced analysis (e.g. Bayesian analysis of the search, estimate, and testing the hypothesis following Gregory & Loredo 1992, hereafter GL) and to explore some effects, such as folding into other possible periods (as well half of the frequency of the most significant peak in the periodogram) and checking of the presence of possible asymmetries in the folded light curves over time and various energy ranges, i.e. peak-to-peak and minimum-to-minimum ratios and phase differences (for details of methodical issues, see Hambaryan et al. 2015). Moreover, in order to solve this challenging problem, a study of the distribution of the spin-phases (resulting from the folding of the photon arrival times into above mentioned various frequencies, see Appendix A) has to be performed.

3.1 X-ray timing analysis

First of all, for each pointed observation (Table 1), we used the cleaned photon arrival times (see Sect. 2) to apply the $\mathrm{Z^{2}_{m}}$ test for periodicity detection in the frequency range of 0.05$-$0.5 Hz in the broad energy band (0.16$-$1.6 keV). Note, that this range includes, as reported in the literature, the spin-frequency 0.119 Hz of RX J0720.4$-$3125 as well as the second harmonic and subharmonic frequencies. As already mentioned in Sect. 3, the periodograms revealed a very significant peak at the frequency 0.119 Hz. Apart from that they also show a less significant second peak at the frequency of 0.0596 Hz, i.e. close to half the frequency of the most significant peak. In Fig. 2, we present the $\mathrm{Z^{2}_{m}}$ statistic periodogram of two observations (Obsids: 0124100101 and 0164560501, see Table 1 and Fig. 1), corresponding to the first and the brightest phase of RX J0720.4$-$3125 in the long-term light curve.

If we apply the H-test developed by de Jager et al. (1989); de Jager & Büsching (2010)333The H-test $H\equiv\mathrm{Z^{2}_{m}-4m+4}$ is a modified $Z_{m}^{2}$ test including higher harmonics and finding optimal number of them to be included in the search and detection of a periodic signal. to the data sets, we identify a second peak in all pointed observations of RX J0720.4$-$3125 also in different energy bands. Moreover, for some energy bands, they do show a more significant peak corresponding to 0.0596 Hz rather than 0.119 Hz.

These results led us to consider the frequency at 0.0596 Hz as the most probable true spin period of RX J0720.4$-$3125 and to perform a detailed analysis of the light curve shapes in different energy bands.

For this purpose, (together with $\mathrm{Z^{2}_{m}}$ and $H$-statistics) we applied the GL method to the different data sets in various energy bands for the frequency ranges 0.05$-$0.07 Hz and 0.10$-$0.14 Hz separately. In general, the GL method first tests if a constant, variable or periodic signal is present in a data set (given photon arrival times of events or a binned data set obeying Poissonian counting statistics) when we have no specific prior information about the shape of the signal.

In the GL method, periodic models are represented by a signal of folded photon arrival times into trial frequency with a light curve shape as a stepwise function with $m$ phase bins per period plus a noise contribution. With such a model we are able to approximate a phase-folded light curve of any shape. Hypotheses for detecting periodic signals represent a class of stepwise, periodic models with the following parameters: trial period, phase, noise parameter, and number of bins ($m$). The most probable model parameters are then estimated by marginalization of the posterior probability over the prior specified range of each parameter. In Bayesian statistics, posterior probability contains a term that penalizes complex models (unless there is no significant evidence to support that hypothesis), so we calculate the posterior probability by marginalizing over a range of models, corresponding to a prior range of number of phase bins, $m$, from 2 to 20. For all cases having strong support for the presence of a periodic signal, using posterior probability density functions we have estimated the most probable values and uncertainties of all parameters, i.e. frequency444In the GL method, a posterior probability density function is evaluated for any parameter of the model by application of Bayes theorem, e.g. for the trial frequency ($\nu$):

$p(\nu|D,M_{m})=\frac{C}{\nu}\int_{0}^{2\pi}d\phi\frac{1}{W_{m}(\nu,\phi)},$

where $C$ and $W_{m}(\nu,\phi)=\frac{N!}{n_{1}!n_{2}!\cdot\cdot\cdot n_{m}!}$ are the normalization constant and number of ways the binned distribution could have arisen “by chance”. $n_{j}(\nu,\phi)$ is the number of events falling into the $j\,\rm{th}$ of $m$ phase bins given the frequency, $\nu$, and phase, $\phi$. $N$ is the total number of photons (for details, see Gregory & Loredo, 1992, 1993). For its uncertainty the highest posterior density (HPD) interval is determined as a probabilistic region around a posterior mode or moment, and is similar to a confidence interval in classical statistics. , phase and the most probable number of phase bins of a periodic signal. Results are presented in Table 2 and Figs. 3 and 4. Using these results, we created phase-folded light curves for all observations of RX J0720.4$-$3125 in three different energy bands (soft: 0.16$-$0.38 keV, hard: 0.38$-$1.6 keV, broad: 0.16$-$1.6 keV, Figs. 10, 11, 12, 13).

In addition, we constructed combined phase-folded light curves (Figs. 5, 14) for all observations and different data groups (see further Sect. 3.2) which clearly show the dependence of the double-peak light curve profile on energy band (cf, e.g Hambaryan et al. 2011). Indeed, the relative strength of the modulation666We have employed the relative strength of modulation

$A=\frac{\sum_{i}|{\rm CR}_{i}-{\rm CR}_{\rm mean}|}{\sum_{i}{\rm CR}_{i}},$

(1)

rather than the pulsed-fraction $PF=({\rm CR}_{\max}-{\rm CR}_{\min})/({\rm CR}_{\max}+{\rm CR}_{\min})$, as an appropriate descriptor of the pulsed emission for complex shape light curves (the strength of variation are independent of the light curve shape, i.e. single or double pulse profile). Here ${\rm CR}$ is the count-rate. Note that our defined relative strength of the modulation for a single-peaked sinusoidal light curve differs by a factor of $2/\pi$ from the definition of semi-amplitude or pulsed-fraction). for all the combined data groups are $0.07\pm 0.02$, $0.05\pm 0.02$, $0.07\pm 0.02$, and $0.12\pm 0.04$ in the energy bands of 0.16$-$0.5 keV, 0.5$-$0.6 keV, 0.6$-$0.7 keV, 0.7$-$2.0 keV, respectively. Thus, the X-ray phase-folded light curves of RX J0720.4$-$3125 (cf, RBS1223 with the spin-period of 10.31 s, Hambaryan et al. 2011) have a markedly double-humped shape with different ratios of maxima in different spectral bands. Moreover, they can be parameterized with a mixed von Mises distribution (see App. A, Fig. 9, Table 5), i.e. with double-peaked light curve profile showing a dependence of the estimated parameters (mean directions, concentrations, and proportion) upon the energy band, indicating that radiation emerges from at least two emitting areas. These spin phase-folded light curves and the estimated parameters from the distribution of phases may serve for rough estimates of the viewing geometry and physical characteristics of the emitting areas (see below).

3.2 Rotational phase-resolved X-ray spectral analysis

In order to create high S/N ratio spin phase-resolved spectra of RX J0720.4$-$3125 we coadded spectra for some sequential pointed observations where it shows approximately similar brightness in the broad energy band (0.16-1.6 keV). Namely, the total observational interval, spanned almost 12.5 years, was divided into five groups (see, Fig. 1 vertical dashed-lines) – the second group is representing the brightest phase of RX J0720.4$-$3125 in the hard energy band.

First, we have simultaneously fitted the phase-averaged, high signal-to-noise ratio ($>$ 30, each spectral bin containing at least 1000 counts) spectral data, collected from different observations (Table 1), using a combination of an absorbed blackbody and a Gaussian absorption line (a multiplicative component) in the model (in XSPEC tbabsbbodyradgabs, see Fig. 6 and Table 3). During the fitting the interstellar absorption and the Gaussian line energy was kept linked or free (see, Table 3), while other parameters (blackbody temperature, Gaussian width, and normalizations) were left free for different data groups. Note, that if the linked parameters also left free, the fitting results remained almost unchanged. It is worthwhile to note that these spectra cannot be fitted with the model of fixed parameters of the absorption line feature at the pre-brightening phase.

Next, we have extracted high signal-to-noise ratio ($>$ 10, each spectral bin containing at least 100 counts) spectral data corresponding to the 8 rotational phase bins (the maxima, minima, and intermediate stages) of the double-peaked pulse profile and collected from different observations (see, Fig. 5 14). These 40 spectra for the five data groups were subject of the simultaneous fitting with a number of highly magnetized INS surface/atmosphere models developed by us (Suleimanov et al. 2010) and implemented into the X-ray spectral fitting package XSPEC (Hambaryan et al. 2011). They are based on various local surface models and compute rotational phase dependent integral emergent spectra of INS, using analytical approximations. The basic model includes temperature/magnetic field distributions over the INS surface, viewing geometry and gravitational redshift. Three local radiating surface models are also considered, namely a naked condensed iron surface and partially ionized hydrogen model atmosphere, semi-infinite or finite on top of the condensed surface. To compute an integral spectrum, the model uses an analytical expression for the local spectra, an emitting condensed iron surface and a diluted blackbody spectrum, transmitted through the thin magnetized atmosphere with an absorption feature. The center of the line depends on the local magnetic field strength, representing either the ion cyclotron line for completely ionized iron and/or the blend of the proton cyclotron and bound-bound atomic transitions in neutral hydrogen. In the latter case, it is parameterized by half of a Gaussian line with parameters: $E_{\mathrm{line}}$, optical depths $\tau=\tau^{0}\exp\left(-\frac{(E-E_{\mathrm{line}})^{2}}{2\sigma_{\mathrm{line}}^{2}}\right)$ and the widths $\sigma_{\mathrm{line}}$ (for details and further references, see Suleimanov et al. (2010); Hambaryan et al. (2011)).

Given the large number of free parameters in the model, we performed a preliminary analysis in order to get rough estimates (or constraints) on some of them. From the observed double-peaked profile light curve shape in different energy bands (see, right panel Fig. 5 and Fig. 14), it has been already clear that the two emitting areas have slightly different spectral and geometrical characteristics (e.g., the relatively cooler one has a larger size). Secondly, some constraints (initial values for the fitting) on magnetic field strengths at the poles may be used on the base of period and its derivative values999A steady spin-down of RX J0720.4$-$3125 was derived (Hambaryan & Neuhäuser 2016) based on the spin-frequency history (Fig. 4) by robust linear fitting, i.e. eliminating the outliers : $\mathrm{P_{0}=16.78215s,\dot{P}=1.86\times 10^{-13}ss^{-1},MJD=51677.104103}$, and $\mathrm{\Delta\nu\sim 10^{-6}-10^{-5}Hz}$ at MJD=53147.442236 assuming magnetic dipole braking as a main mechanism of the spin-down of RX J0720.4$-$3125 (see, e.g. Hohle et al. 2012b; Hambaryan & Neuhäuser 2016). Moreover, from the peak-to-peak separation and the ratio of the minima in this double-humped light curve, we may locate the cooler one at some offset angle with respect to the magnetic axis and azimuth (Schwope et al. 2005).

Having these general constraints on temperatures and magnetic field strengths at the poles, we simulated a large number of photon spectra (absorbed blackbody with Gaussian absorption line) folded with the response of XMM-Newton  EPIC-pn camera with the same observing mode of RX J0720.4$-$3125. We included also the interstellar absorption (for parameters see Table 3), and used a characteristic magnetic field strength of B$\sim$ 3.0$\times$1013 Gauss. The predicted phase-folded light curves (with a purely dipolar magnetic field configuration, depending on the relative orientation between the magnetic axis, the rotation axis, and the observer line of sight, as well gravitational redshift, see upper left panel of Fig. 9 in Appendix A) in four spectral ranges (0.16-0.5 keV, 0.5-0.6 keV, 0.6-0.7 keV, 0.7-2.0 keV) normalized to the maximum of the brightness, were cross-correlated with observed ones (Fig. 5). We infered some constraints on the parameters from the unimodal distributions of them, when cross-correlation coefficients were exceeding 0.9 in the mentioned four energy bands simultaneously and used them as initial input value and as lower and upper bounds for fitting purposes. For example, gravitational redshift cannot exceed 0.3 (it is not possible to obtain the observed relative strength of modulation at larger $z$ due to strong light bending) or the antipodal shift angle must be less than $5^{\circ}$ (due to the observed phase separation between the two peaks in the light curve). It is also evident that the sum of the inclination angle of the line of sight and magnetic poles relative to the rotational axis are already constrained by the light curve class (see, Poutanen & Beloborodov 2006, class III).

Having the above-mentioned crude constraints and input values of the free parameters, we performed fitting with the models implemented in the XSPEC package of all the combined spectra of RX J0720.4$-$3125 (see Table 1 and Fig. 1) simultaneously, i.e. each of those phase resolved spectra considered as different data set with the fixed spin-phase, and with the linked parameters to the other data sets.

The fitted model (Chi$-$Squared fit statistic value of 6070 for 4910 degrees of freedom) with dipolar field configuration (magnetic field strengths at the poles of $\mathrm{B_{p1}\approx B_{p2}\approx(0.3\pm 0.02)\times 10^{14}\leavevmode\nobreak\ G}$), as an orthogonal rotator (angle between the magnetic and the rotation axis $\theta=(89\pm 5)^{\degr}$), with the inclination angle (angle between rotation axis and the observer line of sight) of $i=(42\pm 3)^{\degr}$, and the shift angle of $\kappa\approx 0^{\degr}$ is presented in Table 4 (see also, upper left panel of Fig. 9 in Appendix A).

In order to assess the degree of uniqueness and to estimate confidence intervals of the determined parameters, we have additionally performed a Markov Chain Monte Carlo (MCMC) fitting as implemented in XSPEC. The parameters of the fitted model are presented in Table 4 and Fig. 7.

In Fig. 8 we present the probability density distributions of the gravitational redshift. Note, independent initial input values of parameters of the MCMC approach converged to the same values, in different chains.

4 Discussion

Our comprehensive X-ray timing analysis of RX J0720.4$-$3125, based on the XMM-Newton EPIC-pn long-term pointed observations, showed that the most probable, plausible spin frequency of this INS can be considered 0.0596 Hz (instead of 0.119 Hz reported in the literature). The phase-folded light curves of this period show a markedly double-peaked shape, which formally can be parameterized with a mixed von Mises distribution with significantly different components (more precisely, concentration and mixing proportion parameters, see Appendix A).

A similar picture we had seen with RBS 1223, another isolated neutron star. The pulsation discovery paper reported on a spin-period 5.16 sec (Hambaryan et al. 2002); however, a more detailed study showed that true period is 10.31 s (Haberl et al. 2003; Haberl 2004; Schwope et al. 2005, 2007). Namely, this period with its double peaked profile was suggested by the presence of a second peak (with lower significance) in the power spectrum (Haberl 2004). Further analysis of light curves in different energy bands and phase resolved spectroscopic study (Schwope et al. 2007; Hambaryan et al. 2011) provided additional confirmation on a plausibility of the model of RBS 1223 with emergent emission dominated by two emitting areas with slightly different parameters. Yet another case, Speagle et al. (2011) noted that without the radio period, they might have actually identified the pulse period of PSR J0726 as half of the true value. Cropper et al. (2001) performed spin pulse profile analysis of RX J0720.4$-$3125 based on earlier XMM-Newton  observations, also folded photon arrival times with twice of the detected period via Fourier periodogram which produced a double-humped light curve with different count rates at the different peak; however, having no indication of subharmonics in the periodogram, they did not really consider twice of the period as a true one. Another possible indication of the double-peaked light curve shape of RX J0720.4$-$3125 may be the non-sinusoidal light curve shape and the reported phase lag in different energy bands (Haberl et al. 2006; Hohle et al. 2009). Indeed, slightly different effective temperatures/sizes of two emitting areas at the magnetic poles, located not exactly podal and antipodal directions, may cause such kind of lag, if photon arrival times were to fold into the phases with half of the value of the genuine spin-period (see Fig. 9 in Appendix A).

The combined phase-resolved spectra of RX J0720.4$-$3125 can be simultaneously fitted by emergent radiation of a spectral model of an iron condensed surface with a partially ionized hydrogen atmosphere above it. Formally they can be fitted also by a blackbody spectrum with proton-cyclotron absorption Gaussian line and a peaked (typical for an electron scattering atmosphere) angular distribution of the emergent radiation. In both cases two emitting areas with slightly different characteristics are required (see Table 4).

It is worth to note, that the resulting fit parameters are very similar for different spectral models, also confirmed by an MCMC approach with different input parameters.

However, we believe that the emission properties due to the condensed surface model with a partially ionized, optically thin hydrogen layer above it, including vacuum polarization effects, is more physically motivated. Moreover, semi-infinite atmospheres have rather fan-beamed emergent radiation (see Suleimanov et al. 2010; Hambaryan et al. 2011, and references therein) and it seems impossible to combine a proton cyclotron line with a pencil-beamed emergent radiation.

Emission spectra based on realistic temperature and magnetic field distributions with strongly magnetized hydrogen atmospheres (or other light elements) are formally still an alternative121212Our attempt to fit the combined, phase-averaged spectrum of RX J0720.4$-$3125 by partially ionized, strongly magnetized hydrogen or mid-Z element plasma model (XSPEC nsmax, Mori & Ho 2007; Ho et al. 2008), as well two spots or purely condensed iron surface models, failed. Noteworthy, an acceptable fit is obtained by nsmax model with additional, multiplicative Gaussian absorption line component (model gabs)., but it is unphysical because the absorption line is added by hand.

A pure proton-cyclotron absorption line scenario can be excluded owing to the width of the observed absorption spectral feature in the X-ray spectrum of RX J0720.4$-$3125. Magnetized semi-infinite atmospheres predict too low an equivalent width of the proton cyclotron line in comparison with the observed one.

This result of spectral modelling (i.e the fitting with the condensed surface model with partially ionized, optically thin hydrogen atmosphere above it, including vacuum polarization effects) suggests a true radius of RX J0720.4$-$3125 of $13.3\pm 0.5\mathrm{km}$ for a standard neutron star of 1.4 solar mass, and indicates a stiff equation of state of RX J0720.4$-$3125 (for similar results, see also Ho et al. 2007; Heinke et al. 2006; Suleimanov & Poutanen 2006; Suleimanov et al. 2011).

Our model fits provide the angles between our line of sight and both the rotation axis and the magnetic field axes, we obtain $42\pm 3^{\circ}$ and $47\pm 6^{\circ}$, respectively. The direction of motion should be similar to the former angle in case of spin-orbit almost (see, e.g., Ng & Romani 2007, $\theta<15^{\circ}$) alignment.

In Tetzlaff et al. (2011), the direction of motion was determined to be $62\pm 9^{\circ}$, by tracing back the motion of RX J0720.4$-$3125 to its possible birth association Trumpler$-$10. The direction of motion determined kinematically is thus roughly consistent with almost alignment of spin and orbit.

5 Conclusions

We have carried out a complete X-ray timing re-analysis of the multi-epoch observational data sets of RX J0720.4$-$3125 for an identification of the true spin-period. The genuine spin-period of the isolated neutron star RX J0720.4$-$3125 is twice of the previously reported value in the literature (16.78 s instead of 8.39 s), with a markedly double-peaked light curve profile depending on time and energy.

The observed phase-resolved spectra of the INS RX J0720.4$-$3125 are satisfactorily fitted with slightly different physical and geometrical characteristics of two emitting areas, for the model of a condensed iron surface, with partially ionized, optically thin hydrogen atmosphere above it, including vacuum polarization effects, as orthogonal rotator. The fit also suggests the absence of a strong toroidal magnetic field component. Moreover, the determined mass-radius ratio ($(M/M_{\sun})/(R/\mathrm{km})=0.105\pm 0.002$) suggests a very stiff equation of state of RX J0720.4$-$3125.

More work for detailed spectral model (also with quadrupole magnetic field configuration) computation will be certainly worth to do in the near future and its application to the phase-resolved spectra of other INSs and magnetars. In particular, including also high resolution spectra observed by XMM-Newton and Chandra with possibly other absorption and emission features (Hambaryan et al. 2009; Potekhin 2010; Hambaryan & Neuhäuser 2016).

Acknowledgements.

VH, VS, and RN acknowledge support by the German National Science Foundation (Deutsche Forschungsgemeinschaft, DFG) through project C7 of SFB/TR 7 “Gravitationswellenastronomie”.

Appendix A True spin-period of RX J0720.4$-$3125

Most of the methods for the search and detection of a periodicity in X-ray astronomy are based on the analysis of the distribution of the phases ($\theta_{i}=\nu\,(t_{i}-t_{0})+\dot{\nu}\,(t_{i}-t_{0})^{2}/2$), which is the result of folding photon arrival times ($t_{i}$ ($i=1,...,N$) with a trial frequency ($\nu$) and its derivative ($\dot{\nu}$) for the epoch of zero phase ($t_{0}$, e.g. $\mathrm{Z^{2}_{m}}$ , H-test, the GL method, etc.).

This procedure transforms the time variable to the phases distributed on a circle.

In order to disentangle between single- or multiple-peaked profile of a light curve shape it has to be analyzed statistically. For this purpose, commonly used circular distributions are the von Mises distribution, the wrapped normal distribution, the wrapped Cauchy distribution, etc. The von Mises distribution ($f(\theta;\mu,\kappa)=\frac{1}{I_{0}(\kappa)}e^{\kappa\cos(\theta-\mu)}$) has two parameters: the mean direction and concentration of the data distributed on a circle ($I_{0}(\kappa)$ is the modified Bessel function of order 0).

For demonstration purposes, in Figure 9, we show a sample for the light curve shape modeled by an exponentiated Fourier series, or generalized von Mises distribution (Connors 1997) for two component ($r(t_{i})\propto e^{\sum_{j=1}^{2}{[-\kappa_{j}\cos(\frac{2\pi t_{i}}{P}+\varphi_{j})]}}$, with phases $\varphi_{1,2}$ of photon arrival times ($t_{i},i=1,N$) of a periodic signal with the period of $P$) and distribution of phases (rose diagram-a circular histogram plot of a mixed two von Mises distribution (Mardia & Jupp 2009; Connors 1997), i.e. $f(\theta;\mu_{1},\mu_{2},\kappa_{1},\kappa_{2},p)=\frac{p}{I_{0}(\kappa_{1})}e^{\kappa_{1}\cos(\theta-\mu_{1})}+\frac{1-p}{I_{0}(\kappa_{2})}e^{\kappa_{2}\cos(\theta-\mu_{2})}$, where $\mu_{1},\mu_{2}$ are the mean directions of the phases, $\kappa_{1},\kappa_{2}$ are concentrations, and $p$ is a mixing proportion. These parameters are determining the locations, peakedness, relative contribution of each component of a double-profile light curve, i.e. are expressing the effective temperatures and sizes of the dominating areas resulting the emergent radiation.

To fit the double-peaked light curve profile we coded an Expectation-Maximization algorithm (see, e.g. Mardia et al. 2007, 2008, an iterative approach to find the maximum of the likelihood, applied to a mixture models) used to model rotational phases of all registered photons (column $N$ in the Table 5) as a mixture of two component von Mises distribution. The results of the fit are presented in the Table 5, where for the estimation of uncertainties, an additional, bootstrap algorithm was implemented and applied. We generate a bootstrap sample of 30 datasets for each data group, by randomly sampling from original spin-phases without replacement. Each of this bootstrap datasets fitted with a mixture of two component von Mises distribution. The standard deviations of the derived parameters accepted as an uncertainty measure.

An application of so-called one sample t-test to the difference of the locations of the peaks (column $\mu_{1}/\mu_{2}$, Table 5, or two-sample paired t-test to the $\mu 1$ and $\mu 2$) of the phase-folded light curve (double-peaked profile) is consistent to be about $\pi$, regardless of the energy band. Hence, the null hypothesis on the equality of the ratio of the mean directions of 0.5 (i.e. podal and antipodal directions) is accepted. In contrary, for the parameters of concentration and proportion ($\kappa_{1}/\kappa_{2}$ and $p$) the null hypothesis is rejected ($\kappa_{1}$ and $\kappa_{2}$ are nonidentical populations and $p$ mean value is not equal to 0.5, see, Table 5). We arrive to the similar conclusions by an application of a non-parametric Wilcoxon test.

Appendix B Phase-folded light curves

Here (Fig. 10), we present rotational phase-folded light curves for all pointed observations of RX J0720.4$-$3125 performed by XMM-Newton EPIC-pn in different energy bands (0.16-0.38 keV, 0.38-1.6 keV, and 0.16-1.6 keV) with spin-frequencies (Table 2), zero-phases (derived by an application of the GL method), and using 20 phase bins. All double-peaked profile light curves are normalized to the mean count rate of the pointed observations in the corresponding energy band.

Appendix C Phase-resolved spectra with MCMC

As already mentioned, in order to assess the degree of uniqueness and to estimate confidence intervals of the determined parameters for the model of an emitting condensed iron surface with partially ionized hydrogen model atmosphere, (see, Section 3.2), we have additionally performed a Markov Chain Monte Carlo (MCMC) fitting using the Goodman-Weare chain generating algorithm (the default in XSPEC). The results for the basic parameters of the fitted model (for two data groups I and II) are presented as a matrix plot (see, Fig. 15).