A deep XMM-Newton look on the thermally emitting isolated neutron star RX J1605.3+3249

TL;DR

This study used extensive XMM-Newton observations to search for periodic signals and spectral features in the isolated neutron star RX J1605.3+3249, ultimately setting strict upper limits on pulsations and refining its spectral model.

Contribution

The paper provides the first deep timing upper limits on RX J1605.3+3249, refines its spectral parameters, and investigates the absence of previously reported features, advancing understanding of its emission properties.

Findings

No periodic modulation above 1.5% (3σ) was detected.

Spectral analysis favors a neutron star hydrogen atmosphere with a broad absorption line.

No evidence found for the previously reported narrow absorption feature at 570 eV.

Abstract

Previous XMM-Newton observations of the thermally emitting isolated neutron star RX J1605.3+3249 provided a candidate for a shallow periodic signal and evidence of a fast spin down, which suggested a high dipolar magnetic field and an evolution from a magnetar. We obtained a large programme with XMM-Newton to confirm its candidate timing solution, understand the energy-dependent amplitude of the modulation, and investigate the spectral features of the source. We performed extensive high-resolution and broadband periodicity searches in the new observations, using the combined photons of the three EPIC cameras and allowing for moderate changes of pulsed fraction and the optimal energy range for detection. A deep $4\sigma$ upper limit of $1.33(6)\%$ for modulations in the relevant frequency range conservatively rules out the candidate period previously reported. Blind searches revealed no other periodic signal above the $1.5\%$ level $(3\sigma$) in any of the four new observations. While theoretical models fall short at physically describing the complex energy distribution of the source, best-fit X-ray spectral parameters are obtained for a fully or partially ionized neutron star hydrogen atmosphere model with $B=10^{13}$ G, modified by a broad Gaussian absorption line at $385\pm10$ eV. The deep limits from the timing analysis disfavour equally well-fit double temperature blackbody models where both the star surface and small hotspots contribute to the X-ray flux of the source. We identified a low significance ($1\sigma$) temporal trend on the parameters of the source in the analysis of RGS data dating back to 2002, which may be explained by unaccounted calibration issues and spectral model uncertainties. The new dataset also shows no evidence of the previously reported narrow absorption feature at $\sim570$ eV, whose possible transient nature disfavours an atmospheric origin.

Figures

Tables

Table 1. Table 1: Log of the XMM-Newton AO14 and 2012 observations of RX~J1605.3+3249

obsid	Date	Inst.	Mode	Duration	GTI
				(s)	(%)
0764460201	2015-07-21	pn	FF	118 344	97
	2015-07-21	MOS1	LW	110 089	100
	2015-07-21	MOS2	LW	119 079	100
	2015-07-21	RGS1	SES	120 169	100
	2015-07-21	RGS2	SES	120 113	100
0764460301	2015-07-26	pn	FF	65 038	96
	2015-07-26	MOS1	LW	66 650	100
	2015-07-26	MOS2	LW	66 629	100
	2015-07-26	RGS1	SES	66 858	100
	2015-07-26	RGS2	SES	66 757	100
0764460401	2015-08-20	pn	FF	68 656	81
	2015-08-20	MOS1	LW	70 459	91
	2015-08-20	MOS2	LW	71 629	90
	2015-08-20	RGS1	SES	71 858	91
	2015-08-20	RGS2	SES	71 777	89
0764460501	2016-02-10	pn	FF	59 932	96
	2016-02-10	MOS1	LW	61 540	100
	2016-02-10	MOS2	LW	61 495	100
	2016-02-10	RGS1	SES	61 749	100
	2016-02-10	RGS2	SES	61 675	100
0671620101	2012-03-06	pn	FF	58 542	68
	2012-03-06	MOS1	FF	57 432	89
	2012-03-06	MOS2	FF	57 008	90
	2012-03-06	RGS1	SES	60 406	88
	2012-03-06	RGS2	SES	60 406	84

Table 2. Table 2: Summary of archival XMM-Newton observations of RX~J1605.3+3249  used in the RGS spectral analysis

Ref.	obsid	Date	Net exposure (ks)
			RGS1	RGS2
(A)	0073140201	2002-01-15	27.5	27.0
(B)	0073140301	2002-01-09	18.6	17.0
(C)	0073140501	2002-01-19	22.0	21.5
(D)	0157360401	2003-01-17	29.3	28.5
2012	0671620101	2012-03-06	53.4	50.8
Total archival RGS data (ks)			151	145
Total RGS data⋆⋆footnotemark: ⋆ (ks)			465	457

Table 3. Table 3: Parameters of RX~J1605.3+3249 , as extracted from the EPIC images in the AO14 and 2012 observations

Parameter / obsid	0764460201	0764460301	0764460401	0764460501	0671620101⋆⋆footnotemark: ⋆
Detection likelihood	$2.5\times {10}^6$	$1.4\times {10}^6$	$1.2\times {10}^6$	$1.2\times {10}^6$	$0.9\times {10}^6$
Counts	$4.360\left(7\right)\times {10}^5$	$2.458\left(5\right)\times {10}^5$	$2.185\left(5\right)\times {10}^5$	$2.262\left(5\right)\times {10}^5$	$1.669\left(5\right)\times {10}^5$
…$0.2-0.5$ keV	$2.579\left(5\right)\times {10}^5$	$1.456\left(4\right)\times {10}^5$	$1.305\left(4\right)\times {10}^5$	$1.331\left(4\right)\times {10}^5$	$0.987\left(3\right)\times {10}^5$
…$0.5-1.0$ keV	$1.669\left(4\right)\times {10}^5$	$0.939\left(3\right)\times {10}^5$	$0.823\left(3\right)\times {10}^5$	$0.872\left(3\right)\times {10}^5$	$0.6356\left(27\right)\times {10}^5$
…$1.0-2.0$ keV	$1.11\left(12\right)\times {10}^4$	$0.625\left(9\right)\times {10}^4$	$0.566\left(8\right)\times {10}^4$	$0.586\left(9\right)\times {10}^4$	$0.462\left(8\right)\times {10}^4$
Rate (s-1)	$4.278\left(7\right)$	$4.326\left(9\right)$	$4.226\left(10\right)$	$4.351\left(10\right)$	$4.379\left(12\right)$
Rate (frame-1 camera-1)
…pn	$0.2832\left(4\right)$	$0.2861\left(6\right)$	$0.2804\left(6\right)$	$0.2861\left(6\right)$	$0.2842\left(7\right)$
…MOS1	$0.5601\left(21\right)$	$0.5992\left(28\right)$	$0.5283\left(27\right)$	$0.5921\left(29\right)$	$1.589\left(10\right)$††footnotemark: †
…MOS2	$0.6338\left(22\right)$	$0.6328\left(9\right)$	$0.639\left(3\right)$	$0.641\left(3\right)$	$1.889\left(11\right)$††footnotemark: †
‡‡footnotemark: ‡HR1	$-0.2142\left(15\right)$	$-0.2160\left(21\right)$	$-0.2262\left(23\right)$	$-0.2082\left(22\right)$	$-0.2164\left(26\right)$
‡‡footnotemark: ‡HR2	$-0.8757\left(13\right)$	$-0.8753\left(17\right)$	$-0.8714\left(18\right)$	$-0.8740\left(18\right)$	$-0.8644\left(21\right)$
‡‡footnotemark: ‡HR3	$-0.9975\left(13\right)$	$-0.9958\left(21\right)$	$-0.9993\left(16\right)$	$-0.9994\left(12\right)$	$-0.9955\left(26\right)$
RA (h min sec)	$16$  $05$  $18.5\left(6\right)$	$16$  $05$  $18.5\left(6\right)$	$16$  $05$  $18.5\left(7\right)$	$16$  $05$  $18.5\left(8\right)$	$16$  $05$  $18.4\left(9\right)$
DEC (d m s)	$+32$  $49$  $19.3\left(5\right)$	$+32$  $49$  $19.2\left(5\right)$	$+32$  $49$  $19.7\left(6\right)$	$+32$  $49$  $19.6\left(7\right)$	$+32$  $49$  $18.7\left(8\right)$
RA offset (′′)	$-0.6\pm 0.4$	$+0.2\pm 0.3$	$-0.7\pm 0.4$	$-0.2\pm 0.5$	$-1.2\pm 0.3$
DEC offset (′′)	$+1.1\pm 0.3$	$+1.0\pm 0.3$	$+0.5\pm 0.4$	$+0.9\pm 0.5$	$-0.3\pm 0.5$
Reference sources	$56$	$53$	$50$	$49$	$38$

Table 4. Table 4: Upper limits on pulsations from the timing analysis

	EPIC ($0.01-0.56$ Hz)		pn ( Hz)
Band⋆⋆footnotemark: ⋆	Pulsed fraction
	(%, $3\sigma$ c.l.)
soft	$1.36\pm 0.05$		$1.60\pm 0.10$
hard	$2.01\pm 0.07$		$2.78\pm 0.16$
total	$1.44\pm 0.05$		$1.38\pm 0.08$

Table 5. Table 5: Results of the best-fit double blackbody model with a Gaussian absorption line (per observation and camera)

obsid	${\chi }_{\nu }^2$	NHP	$N_H$a𝑎aa𝑎a$a$The column density is in units of ${10}^{20}$ cm-2.	$k{T}_1^{\infty }$	$k{T}_2^{\infty }$	$R_1^{\infty }$b𝑏bb𝑏b$b$The radiation radius at infinity for each component is computed from the derived blackbody luminosity for a source at a distance of $d\equiv d_{300}=300$ pc.	$R_2^{\infty }$b𝑏bb𝑏b$b$The radiation radius at infinity for each component is computed from the derived blackbody luminosity for a source at a distance of $d\equiv d_{300}=300$ pc.	$ϵ$	$\sigma$	$EW$	$f_X$c𝑐cc𝑐c$c$The observed model flux is in units of ${10}^{-12}$ erg s-1 cm-2 in energy band $0.2-12$ keV.
		(%)		(eV)	(eV)	(km)	(km)	($3\sigma$, eV)	(eV)	($1\sigma$, eV)
pn
2012	$1.18$	$26$	${4.4}_{-3}^{+1.6}$	${64}_{-4}^{+11}$	${119.2}_{-2.4}^{+2.6}$	${13.8}_{-5}^{+1.9}$	${1.22}_{-0.09}^{+0.10}$		$119\left(10\right)$		$6.37\left(5\right)$
201	$0.89$	$88$	${3.4}_{-0.7}^{+1.3}$	${64}_{-4}^{+3}$	${120.0}_{-1.6}^{+1.7}$	${12.2}_{-2.0}^{+0.8}$	$1.15\left(6\right)$		${116}_{-8}^{+4}$		$6.73\left(3\right)$
301	$1.04$	$32$	${5.4}_{-1.0}^{+1.7}$	$59\left(4\right)$	${114.0}_{-1.5}^{+1.7}$	${18.4}_{-3}^{+1.5}$	${1.50}_{-0.07}^{+0.08}$		${114}_{-11}^{+4}$		$6.28\left(3\right)$
401	$0.96$	$63$	${4.2}_{-1.2}^{+1.7}$	${61}_{-5}^{+6}$	${117.6}_{-1.9}^{+2.0}$	${15.7}_{-4}^{+1.6}$	${1.28}_{-0.07}^{+0.08}$		${111}_{-9}^{+5}$		$6.58\left(5\right)$
501	$0.88$	$61$	${5.4}_{-1.3}^{+2.2}$	$57\left(4\right)$	$114.0\left(1.6\right)$	${20.4}_{-3}^{+2.1}$	$1.55\left(6\right)$		${109}_{-12}^{+5}$		$6.46\left(4\right)$
MOSd𝑑dd𝑑d$d$Simple fit of combined MOS1 and MOS2 spectra (per observation).
2012	$1.02$	$41$	${1.5}_{-1.4}^{+1.6}$	${74}_{-10}^{+9}$	$127\left(3\right)$	${6.3}_{-2.4}^{+0.8}$	${0.84}_{-0.06}^{+0.07}$		${95}_{-16}^{+20}$		$6.58\left(6\right)$
201	$1.09$	$19$	${1.8}_{-0.7}^{+1.2}$	${70}_{-5}^{+4}$	$125.6\left(2.2\right)$	${8.1}_{-1.6}^{+0.7}$	$0.88\left(5\right)$		${119}_{-8}^{+3}$		$6.70\left(4\right)$
301	$0.95$	$66$	${2.0}_{-1.4}^{+1.6}$	${77}_{-8}^{+9}$	$131\left(4\right)$	${5.8}_{-1.9}^{+0.5}$	${0.73}_{-0.06}^{+0.07}$		${91}_{-15}^{+17}$		$6.29\left(4\right)$
401	$1.05$	$31$	${2.1}_{-1.0}^{+1.4}$	$69\left(5\right)$	$126\left(3\right)$	${9.0}_{-1.6}^{+0.9}$	$0.85\left(6\right)$		${120}_{-13}^{+4}$		$6.58\left(5\right)$
501	$1.06$	$22$	${3.6}_{-2.3}^{+2.2}$	${43}_{-5}^{+8}$	${118.3}_{-1.6}^{+1.7}$	$>20$	$1.26\left(6\right)$	$520\left(30\right)$	${121}_{-12}^{+13}$		$7.01\left(6\right)$
EPICe𝑒ee𝑒e$e$Simultaneous fit of pn and combined MOS spectra (per observation).
2012	$1.15$	$2.2$	${3.9}_{-0.9}^{+1.0}$	${64.2}_{-2.8}^{+6}$	$120.5\left(1.8\right)$	${13.3}_{-2.9}^{+1.0}$	$1.16\left(5\right)$		${117}_{-16}^{+3}$		$6.40\left(3\right)$
201	$1.00$	$48$	${2.4}_{-0.9}^{+0.5}$	${67.3}_{-1.7}^{+3}$	$122.3\left(1.2\right)$	${10.4}_{-1.2}^{+0.5}$	$1.07\left(3\right)$		${118.6}_{-4}^{+2.5}$		$6.885\left(25\right)$
301	$1.02$	$37$	${3.2}_{-1.3}^{+1.4}$	${68}_{-5}^{+6}$	${120.7}_{-1.6}^{+1.7}$	${9.5}_{-2.3}^{+0.6}$	$1.12\left(5\right)$		${103}_{-13}^{+15}$		$6.45\left(3\right)$
401	$1.01$	$41$	${2.9}_{-0.7}^{+1.1}$	${65}_{-3}^{+4}$	$121.6\left(1.8\right)$	${11.3}_{-1.9}^{+0.8}$	$1.06\left(5\right)$		${116}_{-10}^{+3}$		$6.77\left(3\right)$
501	$1.08$	$14$	${3.0}_{-0.6}^{+1.2}$	${66.4}_{-6}^{+1.2}$	$123.6\left(1.8\right)$	${10.9}_{-2.1}^{+0.7}$	${1.01}_{-0.04}^{+0.05}$		${118}_{-7}^{+3}$		$6.82\left(3\right)$
Multi-epoch fitsf𝑓ff𝑓f$f$Simultaneous fit per instrument (5 pn, 10 MOS1/2, and 6 pn and stacked MOS spectra).
pn	$1.02$	$29$	${4.5}_{-0.6}^{+0.5}$	${60.9}_{-1.5}^{+1.7}$	$117.0\left(8\right)$	${16.2}_{-1.3}^{+0.6}$	$1.34\left(3\right)$		${114.6}_{-2.6}^{+2.0}$	${100}_{-60}^{+50}$	$6.468\left(16\right)$
MOS	$1.11$		${2.5}_{-0.4}^{+0.7}$	${68.4}_{-3}^{+2.0}$	$125.7\left(1.2\right)$	${9.7}_{-1.1}^{+0.4}$	$0.90\left(3\right)$		${118.9}_{-5}^{+1.6}$	${100}_{-40}^{+50}$	$6.877\left(22\right)$
EPIC	$1.12$		${3.1}_{-0.3}^{+0.5}$	${65.8}_{-2.2}^{+1.0}$	$121.7\left(7\right)$	${11.8}_{-0.9}^{+0.3}$	$1.101\left(20\right)$		${118.1}_{-2.8}^{+1.3}$	${90}_{-10}^{+40}$	$6.663\left(12\right)$
Multi-epoch fitsg𝑔gg𝑔g$g$Results of simultaneous fits of a single temperature model, tbabs(bbody-gauss).				$k{T}^{\infty }$ (eV)		$R^{\infty }$ (km)
pn	$3.8$	$\ll 1$	${2.4}^{\star }$	${88.16}_{-0.19}^{+0.18}$		${5.16}_{-0.04}^{+0.03}$			$170.7\left(8\right)$	${135}_{-19}^{+6}$	$6.077\left(7\right)$
MOS	$2.5$	$\ll 1$	${2.4}^{\star }$	${91.12}_{-0.29}^{+0.3}$		${4.39}_{-0.05}^{+0.04}$			$180.8\left(1.5\right)$	${117}_{-14}^{+26}$	$5.678\left(10\right)$
EPIC	$5.5$	$\ll 1$	${2.4}^{\star }$	${89.18}_{-0.16}^{+0.14}$		$4.94\pm 0.03$			$174.1\left(6\right)$	${112.9}_{-1.1}^{+19}$	$5.928\left(5\right)$

Table 6. Table 6: Results of the best-fit fully ionized neutron star atmosphere model with a Gaussian absorption line (per observation and camera)

obsid	${\chi }_{\nu }^2$	d.o.f	NHP	$N_H$	$T_{\mathrm{eff}}$	$d$	$ϵ$	$\sigma$	$EW$	$f_X$a𝑎aa𝑎a$a$The observed model flux is in units of ${10}^{-12}$ erg s-1 cm-2 in energy band $0.2-12$ keV.
			(%)	(${10}^{20}$ cm-2)	(${10}^5$ K)	(pc)	(eV)	(eV)	(eV)
pn
2012	$1.19$	184	$4$	$3.0\left(4\right)$	${5.69}_{-0.08}^{+0.09}$	$116\left(8\right)$	${380}_{-30}^{+23}$	${95}_{-11}^{+14}$	${70}_{-40}^{+3}$	$6.622\left(23\right)$
201	$0.90$	204	$84$	${2.22}_{-0.26}^{+0.25}$	$5.83\left(5\right)$	$130\left(5\right)$	${376}_{-18}^{+15}$	${90}_{-7}^{+8}$	${55}_{-4}^{+5}$	$6.965\left(16\right)$
301	$1.07$	195	$25$	$2.9\left(3\right)$	$5.68\left(7\right)$	$116\left(6\right)$	${388}_{-22}^{+17}$	${88}_{-9}^{+10}$	${66}_{-40}^{+1.0}$	$6.650\left(18\right)$
401	$0.97$	195	$62$	$2.2\left(3\right)$	$5.85\left(7\right)$	${134}_{-7}^{+6}$	${392}_{-23}^{+17}$	${79}_{-9}^{+11}$	${70}_{-30}^{+1}$	$6.912\left(21\right)$
501	$1.01$	196	$45$	${2.37}_{-0.3}^{+0.29}$	$5.92\left(7\right)$	$137\left(7\right)$	${399}_{-19}^{+15}$	${80}_{-10}^{+8}$	${50}_{-6}^{+4}$	$6.909\left(21\right)$
MOSb𝑏bb𝑏b$b$Simple fit of combined MOS1 and MOS2 spectra (per observation).
2012	$1.05$	181	$30$	$2.2\left(5\right)$	${5.82}_{-0.07}^{+0.08}$	${134}_{-9}^{+8}$	${380}_{-25}^{+19}$	${86}_{-9}^{+10}$	${76}_{-30}^{+1.0}$	$6.42\left(3\right)$
201	$1.11$	199	$14$	$1.81\left(4\right)$	${5.79}_{-0.06}^{+0.05}$	${131}_{-6}^{+5}$	${338}_{-25}^{+20}$	${102}_{-8}^{+10}$	${65}_{-25}^{+4}$	${6.741}_{-0.020}^{+0.022}$
301	$0.98$	187	$59$	$3.3\left(4\right)$	$5.70\left(7\right)$	${120}_{-7}^{+6}$	${401}_{-18}^{+14}$	${83}_{-7}^{+9}$	${75}_{-30}^{+2.0}$	$6.081\left(25\right)$
401	$1.08$	185	$23$	$1.9\left(5\right)$	$5.72\left(9\right)$	${128}_{-9}^{+7}$	${353}_{-28}^{+21}$	${97}_{-9}^{+11}$	${75}_{-40}^{+1.0}$	$6.682\left(29\right)$
501	$1.12$	188	$12$	$1.8\left(6\right)$	${5.72}_{-0.09}^{+0.07}$	${124}_{-8}^{+7}$	${335}_{-30}^{+24}$	${102}_{-9}^{+12}$	${76}_{-40}^{+1.0}$	$7.02\left(3\right)$
EPICc𝑐cc𝑐c$c$Simultaneous fit of pn and combined MOS spectra (per observation).
2012	$1.17$	370	$1.3$	${2.63}_{-0.3}^{+0.29}$	$5.78\left(6\right)$	$124\left(5\right)$	${386}_{-18}^{+15}$	${87}_{-7}^{+8}$	${55}_{-7}^{+3}$	$6.546\left(18\right)$
201	$1.02$	408	$36$	${2.00}_{-0.21}^{+0.22}$	$5.81\left(4\right)$	${129}_{-4}^{+3}$	${357}_{-15}^{+12}$	${97}_{-5}^{+6}$	${60}_{-3}^{+4}$	$6.976\left(11\right)$
301	$1.03$	387	$31$	$3.03\left(23\right)$	$5.71\left(5\right)$	$129\left(6\right)$	${395}_{-14}^{+11}$	${85}_{-6}^{+7}$	${55}_{-2.0}^{+4}$	$6.471\left(19\right)$
401	$1.03$	385	$34$	$2.03\left(28\right)$	$5.79\left(5\right)$	$129\left(5\right)$	${372}_{-18}^{+14}$	${89}_{-7}^{+8}$	${60}_{-3}^{+6}$	$6.893\left(22\right)$
501	$1.10$	389	$8$	$2.18\left(28\right)$	$5.83\left(5\right)$	${128}_{-6}^{+4}$	${369}_{-18}^{+14}$	${92}_{-6}^{+8}$	${60}_{-3}^{+7}$	$6.974\left(16\right)$
Multi-epoch fitsd𝑑dd𝑑d$d$Simultaneous fit per instrument (5 pn, 10 MOS1/2, and 6 pn and stacked MOS spectra).
pn	$1.04$	994	$17$	$2.45\left(14\right)$	$5.81\left(3\right)$	${128.2}_{-3}^{+2.6}$	${385}_{-9}^{+8}$	$87\left(4\right)$	${54}_{-2.0}^{+5}$	$6.842\left(9\right)$
MOS	$1.13$	1738		${2.13}_{-0.19}^{+0.21}$	${5.763}_{-0.04}^{+0.026}$	${126.7}_{-4}^{+2.6}$	${361}_{-11}^{+9}$	$93\left(4\right)$	${64}_{-2.0}^{+4}$	$7.016\left(13\right)$
EPIC	$1.15$	1202		${2.30}_{-0.11}^{+0.12}$	${5.798}_{-0.026}^{+0.018}$	${127.3}_{-2.4}^{+1.8}$	${373}_{-7}^{+6}$	${91.2}_{-2.8}^{+3}$	${57.7}_{-1.4}^{+2.1}$	$6.805\left(6\right)$

Table 7. Table 7: Results of the best-fit partially ionized neutron star atmosphere model with a Gaussian absorption line (per observation and camera)

obsid	${\chi }_{\nu }^2$	d.o.f.	NHP	$N_H$a𝑎aa𝑎a$a$The observed model flux is in units of ${10}^{-12}$ erg s-1 cm-2 in energy band $0.2-12$ keV.	$T_{\mathrm{eff}}$	$d$	$R_{\mathrm{em}}$	$ϵ$	$\sigma$	$EW$	$f_X$a𝑎aa𝑎a$a$The observed model flux is in units of ${10}^{-12}$ erg s-1 cm-2 in energy band $0.2-12$ keV.
			(%)	(${10}^{20}$ cm-2)	(${10}^5$ K)	($3\sigma$, pc)	(km)	(eV)	(eV)	(eV)
pn
2012	$1.18$	185	$5$	$4.07\left(17\right)$	${6.17}_{-0.06}^{+0.11}$		${9.05}_{-0.07}^{+0.4}$	${382}_{-16}^{+9}$	$89\left(5\right)$	${52}_{-12}^{+3}$	$6.581\left(23\right)$
201	$0.88$	205	$89$	${4.12}_{-0.21}^{+0.25}$	${5.85}_{-0.09}^{+0.07}$		${10.3}_{-0.5}^{+0.8}$	${362}_{-17}^{+12}$	${92}_{-5}^{+7}$	${61.7}_{-1.5}^{+3}$	${6.708}_{-0.012}^{+0.014}$
301	$1.07$	196	$25$	${4.06}_{-0.11}^{+0.10}$	${6.16}_{-0.08}^{+0.04}$		${9.0}_{-0.7}^{+0.5}$	$390\left(6\right)$	$84\left(4\right)$	${50}_{-10}^{+1.0}$	$6.589\left(18\right)$
401	$0.96$	196	$65$	${3.32}_{-0.22}^{+0.4}$	${6.33}_{-0.14}^{+0.3}$		${8.7}_{-0.8}^{+0.6}$	${387}_{-22}^{+11}$	${80}_{-7}^{+11}$	$50\left(5\right)$	${6.828}_{-0.20}^{+0.19}$
501	$1.00$	197	$50$	${3.58}_{-0.21}^{+0.4}$	${6.34}_{-0.14}^{+0.20}$	$113\left(6\right)$	${9.0}_{-0.7}^{+0.6}$	${394}_{-20}^{+11}$	${81}_{-6}^{+10}$	${48}_{-6}^{+2}$	${6.409}_{-0.018}^{+0.016}$
MOSb𝑏bb𝑏b$b$Simple fit of combined MOS1 and MOS2 spectra (per observation).
2012	$1.06$	182	$29$	$3.4\left(4\right)$	${6.30}_{-0.15}^{+0.5}$		${8.7}_{-1.2}^{+0.6}$	${373}_{-17}^{+15}$	$87\left(8\right)$	${63}_{-5}^{+9}$	${6.311}_{-0.03}^{+0.029}$
201	$1.10$	200	$16$	$3.12\left(29\right)$	${6.16}_{-0.14}^{+0.13}$		${9.0}_{-0.7}^{+0.9}$	${325}_{-19}^{+15}$	$104\left(6\right)$	${68}_{-4}^{+3}$	${6.608}_{-0.020}^{+0.021}$
301	$0.98$	188	$54$	${4.30}_{-0.17}^{+0.13}$	${6.21}_{-0.08}^{+0.12}$		${8.82}_{-0.7}^{+0.28}$	$396\left(7\right)$	$82\left(4\right)$	${56}_{-3}^{+4}$	$6.072\left(25\right)$
401	$1.07$	186	$25$	${3.09}_{-0.3}^{+0.21}$	${6.22}_{-0.18}^{+0.17}$		${8.6}_{-0.8}^{+0.7}$	${344}_{-13}^{+16}$	${99}_{-7}^{+6}$	${66}_{-9}^{+7}$	${6.555}_{-0.028}^{+0.027}$
501	$1.12$	189	$12$	${2.84}_{-0.27}^{+0.20}$	${6.25}_{-0.16}^{+0.17}$		${8.6}_{-0.9}^{+0.6}$	${322}_{-12}^{+14}$	${105}_{-6}^{+5}$	${83}_{-12}^{+1.0}$	$6.602\left(27\right)$
EPICc𝑐cc𝑐c$c$Simultaneous fit of pn and combined MOS spectra (per observation).
2012	$1.16$	371	$1.8$	${4.11}_{-0.18}^{+0.12}$	${6.02}_{-0.08}^{+0.04}$		${9.5}_{-0.5}^{+0.4}$	${375}_{-7}^{+10}$	${89}_{-5}^{+4}$	${58}_{-8}^{+3}$	$6.419\left(18\right)$
201	$1.01$	409	$44$	$3.50\left(16\right)$	${6.06}_{-0.09}^{+0.06}$		${9.38}_{-0.28}^{+0.4}$	${343}_{-11}^{+10}$	${99}_{-4}^{+5}$	${63}_{-2.2}^{+3}$	$6.391\left(10\right)$
301	$1.03$	388	$31$	$4.11\left(10\right)$	$6.18\left(7\right)$		${8.95}_{-0.4}^{+0.19}$	$392\left(4\right)$	$83.2\left(2.5\right)$	${49}_{-6}^{+1.0}$	${6.430}_{-0.013}^{+0.015}$
401	$1.02$	386	$39$	${3.23}_{-0.23}^{+0.19}$	${6.27}_{-0.14}^{+0.12}$		$8.7\left(6\right)$	${364}_{-12}^{+11}$	$91\left(6\right)$	$58\left(4\right)$	$6.769\left(16\right)$
501	$1.10$	390	$10$	${3.55}_{-0.22}^{+0.18}$	${6.17}_{-0.13}^{+0.08}$		${9.2}_{-0.6}^{+0.5}$	${357}_{-13}^{+12}$	${95}_{-5}^{+6}$	$62\left(3\right)$	${6.424}_{-0.015}^{+0.013}$
Multi-epoch fitsd𝑑dd𝑑d$d$Simultaneous fit per instrument (5 pn, 10 MOS1/2, and 6 pn and stacked MOS spectra).
pn	$1.03$	995	$26$	$4.85\left(14\right)$	${5.56}_{-0.04}^{+0.03}$	${106}_{-4}^{+6}$	${11.9}_{-0.4}^{+0.6}$	${362}_{-9}^{+8}$	${92}_{-3}^{+4}$	${62.5}_{-1.5}^{+2.5}$	${6.519}_{-0.007}^{+0.008}$
MOS	$1.13$	1739		${3.31}_{-0.13}^{+0.07}$	${6.24}_{-0.08}^{+0.05}$		${8.66}_{-0.08}^{+0.3}$	${353}_{-4}^{+6}$	${95.0}_{-2.8}^{+1.9}$	${64.1}_{-5}^{+1.0}$	${6.875}_{-0.011}^{+0.013}$
EPIC	$1.14$	1203		$3.82\left(10\right)$	${6.04}_{-0.05}^{+0.03}$		${9.41}_{-0.23}^{+0.25}$	${361}_{-4}^{+5}$	${93.6}_{-2.1}^{+2.0}$	${60.5}_{-2.8}^{+0.6}$	${6.227}_{-0.006}^{+0.007}$

Table 8. Table 8: Results of the RGS spectral analysis

obsida𝑎aa𝑎a$a$The observations are labeled as in Tables 1 and 2.	$k{T}_{\infty }$	$R_{\infty }$b𝑏bb𝑏b$b$The radiation radius at infinity is computed for a source distance of $d\equiv d_{300}=300$ pc.	$ϵ$	$EW$	$F_X$c𝑐cc𝑐c$c$The unabsorbed model flux is in units of ${10}^{-12}$ erg s-1 cm-2 in energy band $0.2-12$ keV.
	(eV)	(km)	(eV)
(A)	$86\left(5\right)$	${3.2}_{-1.5}^{+1.8}$	${430}_{-24}^{+30}$	${60}_{-16}^{+20}$	${6.5}_{-1.4}^{+2.2}$
(B)	${85}_{-6}^{+8}$	${3.2}_{-1.9}^{+2.1}$	${400}_{-28}^{+40}$	${45}_{-26}^{+30}$	${6.3}_{-2.2}^{+2.8}$
(C)	${83}_{-6}^{+7}$	${3.5}_{-1.9}^{+2.2}$	${390}_{-24}^{+30}$	${55}_{-28}^{+24}$	${6.8}_{-2.1}^{+2.7}$
(D)	${81}_{-4}^{+5}$	${3.9}_{-1.9}^{+2.1}$	${400}_{-16}^{+20}$	${85}_{-19}^{+5}$	${7.8}_{-1.9}^{+2.3}$
2012	${90}_{-6}^{+7}$	${2.6}_{-1.3}^{+1.5}$	${410}_{-25}^{+35}$	${45}_{-25}^{+29}$	${5.4}_{-1.3}^{+1.7}$
201	$92\left(3\right)$	${2.5}_{-0.9}^{+1.0}$	${430}_{-23}^{+30}$	${45}_{-26}^{+6}$	${5.2}_{-0.7}^{+0.8}$
301	${91}_{-5}^{+6}$	${2.5}_{-1.1}^{+1.2}$	${410}_{-26}^{+40}$	${55}_{-40}^{+6}$	${5.1}_{-1.0}^{+1.3}$
401	${92}_{-5}^{+6}$	${2.4}_{-1.1}^{+1.2}$	${420}_{-30}^{+60}$	${19}_{-15}^{+35}$	${5.0}_{-1.0}^{+1.3}$
501	${96}_{-6}^{+7}$	${2.2}_{-1.0}^{+1.1}$	${410}_{-60}^{+70}$		${4.8}_{-1.0}^{+1.3}$
old	${85.7}_{-2.8}^{+3}$	${3.20}_{-1.2}^{+1.3}$	${416}_{-12}^{+15}$	${50}_{-17}^{+2}$	${6.6}_{-0.9}^{+1.0}$
new	${91.7}_{-2.2}^{+2.4}$	${2.42}_{-0.7}^{+0.8}$	${407}_{-12}^{+14}$	${60}_{-7}^{+12}$	$4.9\left(5\right)$

Table 9. Table 9: Investigation of a narrow absorption feature in RGS1 data

obsida𝑎aa𝑎a$a$The observations are labeled as in Tables 1 and 2.	$S/N$	$k{T}_{\infty }$	$E{W}_1$	${ϵ}_2$	${\sigma }_2$	$E{W}_2$
		(eV)
(A)	45	${84}_{-4}^{+6}$	${60}_{-28}^{+22}$	${571.9}_{-1.7}^{+1.4}$		$8$
(B)	40	${83}_{-4}^{+6}$	$70$	$-$	$-$
(C)	40	${83}_{-4}^{+6}$	$70$	$-$	$-$
(D)	50	${84}_{-4}^{+5}$		$-$	$-$
2012	50	${91.1}_{-6}^{+2.7}$	${50}_{-40}^{+20}$	${552.0}_{-1.4}^{+2.1}$		$2.6$
201	90	${94}_{-4}^{+7}$		$-$	$-$
301	65	${86}_{-3}^{+5}$	$50$	$-$	$-$
401	60	${92}_{-5}^{+8}$				$2.4$
501	60	${96}_{-10}^{+5}$	$30$	$-$	$-$
old	90	${83.9}_{-2.3}^{+2.7}$	$70$	$575\left(4\right)$	${5.1}_{-2.2}^{+4}$	$3$
new	150	$93\left(3\right)$	$35$		${16}_{-10}^{+12}$	$4$

Equations

$$p_{\mathrm{f}}=\frac{C{R}_{\max }-C{R}_{\min }}{C{R}_{\max }+C{R}_{\min }}$$

Full text

††offprints: A. M. Pires

1 Introduction

In the usual scenario of magnetic dipole braking in vacuum, the observed secular lengthening of a pulsar spin period, typically of 3 s every $10^{8}$ yr, is the consequence of the torque exerted by the magnetic field on the rotating neutron star. Despite its simplicity, the model provides useful estimates of the evolutionary state of a pulsar, namely, its available rotational power, characteristic age, and surface dipolar magnetic field strength (Ostriker & Gunn, 1969). For the isolated neutron stars (INSs) in our Galaxy, dipolar field estimates span across five orders of magnitude (Manchester et al., 2005). Among the sources with the highest values are those known as magnetars (see Turolla et al., 2015; Kaspi & Beloborodov, 2017; Coti Zelati et al., 2018; Esposito et al., 2018, for recent reviews).

Magnetars are usually observed through their violent bursts of high energy and slow down at a much faster rate than normal pulsars (up to 1 s every 60 years, for the extreme case of the soft-gamma repeater SGR~1806-20). According to the most favoured interpretation (Thompson & Duncan, 1995, 1996), their complex phenomenology, transient behaviour, and bright quiescent X-ray luminosity, much in excess of that from spin down, can be explained by crustal and magnetospheric effects provoked by the decay and rearranging of the enormous stellar magnetic field (Goldreich & Reisenegger, 1992; Pons et al., 2009; Beloborodov & Li, 2016). As a result of field dissipation and braking, it is expected that an evolved magnetar ($\gtrsim 10^{5}$ yr) will be less active than young ones, and have a longer spin period and higher surface temperature than ordinary pulsars of similar age (e.g. Perna & Pons, 2011; Turolla et al., 2011).

The group of X-ray thermally emitting INSs discovered by ROSAT and dubbed the magnificent seven (M7, see Haberl, 2007; Kaplan, 2008; Turolla, 2009, for reviews) may have evolved from such a channel of pulsar evolution (e.g. Heyl & Kulkarni, 1998; Kaplan & van Kerkwijk, 2009; Popov et al., 2010; Viganò et al., 2013). They consist of a rather unique local group of middle-aged ($\sim 10^{5}-10^{6}$ yr), cooling neutron stars, displaying similar low blackbody temperatures ($kT\sim 45-100$ eV), long spin periods ($P\sim 3-17$ s), and moderately strong dipolar magnetic field strengths ($B_{\rm dip}\sim{\rm few}\times 10^{13}$ G). Unlike other X-ray pulsars, their emission is purely thermal with no signs of magnetospheric activity and is believed to originate directly from the neutron star surface.

The source RX J1605.3+3249, as the third brightest among the M7 INSs (Motch et al., 1999), was consequently visited by XMM-Newton in several occasions during the early years of its science operations (see van Kerkwijk et al., 2004; Haberl, 2007, for details on the past investigations of the source). However, these early observations were not deep enough to allow the detection of the neutron star spin signal. Ensuing a visibility gap of six years, a 60 ks observation performed in 2012 finally provided a candidate spin period for the INS (Pires et al., 2014). The amplitude of the shallow and strongly energy-dependent periodic signal was detected close to the sensitivity limit of the data; only the harder portion of the source spectrum (roughly, 30% of all source events at energies above $0.5$ keV) was found to show a significant modulation at the $4\sigma$ level. Nonetheless, a joint timing analysis around the detected signal at $P\sim 3.39$ s, connecting the 2012 dataset with early XMM-Newton observations of the source, hinted at an unprecedentedly high value of spin down. The inferred dipolar magnetic field of $B_{\rm dip}\sim 7.4\times 10^{13}$ G overlaps the magnetar range and, if confirmed, could rank the highest in the group.

Besides, the analysis of the then available XMM-Newton EPIC data on the source confirmed the evidence of a complex, multi-temperature, energy distribution and the presence of absorption features. In Pires et al. (2014) we described the spectrum of the source using a two-component blackbody model of temperatures $60$ eV and $110$ eV, superposed by one or two Gaussian absorption features at around energies $400$ eV and $860$ eV; best fits were found assuming the Galactic column density in the direction of the source, $2.4\times 10^{20}$ cm*-2* (Kalberla et al., 2005, see Section 3.2, for details). These results motivated us to investigate the INS further.

We were granted a large programme of observation with XMM-Newton (programme ID: 76446) for a total duration of 310 ks and four satellite visits in the AO14 observing cycle. We report here the results of this observational campaign. The paper is organised as follows: in Section 2 we describe the XMM-Newton observations and the data reduction. Analysis and results are presented in Section 3. The implications of our results are discussed in Section 4, with particular emphasis on the properties and recent work on the group of M7 INSs. Conclusions and the summary of results are in Section 5.

2 Observations and data reduction

The XMM-Newton observatory (Jansen et al., 2001) targeted the INS RX~J1605.3+3249 (hereafter J1605) in four occasions between July 2015 and February 2016, using EPIC as the prior instrument for the investigation. Table 1 contains information on the scientific exposures and instrumental configuration of the EPIC-pn (Strüder et al., 2001), EPIC-MOS (Turner et al., 2001), and RGS (den Herder et al., 2001) detectors.

We included in the timing and spectral analysis (Section 3) the past 2012 observation of the source (obsid 0671620101; Table 1). For the spectral analysis of RGS data (Section 3.2.2), we additionally included archival XMM-Newton observations of the source dating back to 2002 that were not severely affected by background flares (see Table 2). We applied as criterion a minimum net exposure of 10 ks to include observations in the RGS analysis. The archival EPIC observations of J1605 performed before the visibility gap (analysed and discussed in Pires et al. 2014 and references therein) were not included due to the high background level and the heterogeneous science operating modes and optical blocking filters that were then adopted (see also Section 3.2.2). All archival observations were processed and analysed consistently with the data from the large programme.

The main trigger behind our programme was to confirm the candidate spin period and, by means of a precise timing solution, measure the neutron star’s spin down rate. This is better achieved through a well-sampled ephemeris; incoherent methods result into much less accurate spin down determinations, while a two-dimensional periodicity search (like the one performed in Pires et al. 2014) suffers from a large parameter space and number of independent trials that have to be covered in the case of observations that are years apart. As a result, the pulsar best ($P,\dot{P}$) solution is determined at a low confidence level.

Therefore, the time intervals between the four observations in AO14 (of $5$, $25\pm 3$, and $175\pm 8$ days) were carefully chosen to coherently connect them in phase with three past pn observations of the source, for a total time span of $\sim 5160$ days (Section 3.1). Assuming the properties of the periodic signal as detected in the 2012 observation (Pires et al. 2014 and Section 3.1), we considered for the feasibility a total count rate of 0.8 s*-1* in the three EPIC cameras for source photons with energy above $0.5$ keV; the individual exposure times were required to ensure a significant detection of the 2012 signal while allowing for moderate changes of pulsed fraction, $p_{\rm f}=2.5\%-5\%$ (within $\pm 1\sigma$). Likewise, we chose to operate the MOS and pn cameras in large-window (LW) and full-frame (FF) imaging modes, to provide sufficient time resolution (0.9 s and 73.4 ms, respectively) to measure the 2012 modulation. We adopted the thin filter for both instruments, given its better response at soft X-ray energies.

We performed standard data reduction with SAS 15 (xmmsas_20160201_1833-15.0.0) using up-to-date calibration files and following the analysis guidelines of each instrument444http://xmm2.esac.esa.int/docs/documents. We processed the EPIC exposures using the SAS meta tasks epchain and emchain and applied default corrections. For RGS, we used the SAS routine rgsproc to process the raw data files and create masks for the source and background regions.

Background flares were registered occasionally during the AO14 observations, usually lasting for less than a few kiloseconds. The percentages of good-time intervals (GTIs), filtering out periods of high background activity, are shown in Table 1 for each scientific exposure and observation. Standard count rate thresholds were adopted for pn and MOS; for RGS, we used the background count rate on CCD 9 and applied a threshold of $0.1$ s*-1* to filter the event lists with the SAS task rgsfilter. On average, data loss is small: of 2% for MOS and RGS and 7% for pn. The observation the worst affected by flares was that performed in August 2015, with a percentage of data loss between 9% and 19% depending on the camera. The total net exposures per camera are 303 ks (MOS1), 312 ks (MOS2), 314 ks (RGS1), 312 ks (RGS2), and 290 ks (pn).

For the analysis of EPIC data, we filtered the event lists to exclude ‘bad’ CCD pixels and columns, as well as to retain the pre-defined photon patterns with the highest quality energy calibration – that are, single and double events for pn (pattern $\leq 4$) and single, double, triple, and quadruple for MOS (pattern $\leq 12$). The source centroid and optimal extraction region, with typical sizes of 140”, were defined with the SAS task eregionanalyse in the $0.3-1.35$ keV energy band for each EPIC camera and observation. Background circular regions of size $60^{\prime\prime}$ to $100^{\prime\prime}$ were defined away from the source, on the same CCD of the target whenever possible.

The detected source count rates, hardness ratios, the pile-up-relevant count rates per frame for each camera and observation, and other parameters based on a maximum likelihood fitting are listed in Table 3, with nominal $1\sigma$ statistical uncertainties. The parameters are determined with the SAS task emldetect on images created for each camera, observation, and energy band (only the combined EPIC results are shown; the X-ray emission of the source is compatible with the background level at energies above 2 keV). For comparison, we also list the source parameters as determined from the 2012 EPIC exposure.

Following the guidelines of Jethwa et al. (2015), we ensured that the pile-up levels at aimpoint were within tolerant guidelines for both EPIC detectors. Based on the source spectrum and on the number of counts per frame in each camera (listed in Table 3 for direct comparison with Fig. 5 of Jethwa et al. 2015), we estimate that the percentage levels of spectral distortion and flux loss were around 1.5% and 3.5% for pn, and 0.3% and 1% for MOS.

Overall, the source properties are consistent between epochs since 2012: potential discrepancies can be asserted to the cross-calibration uncertainties between the EPIC instruments and to different background levels. In Section 3.2 we investigate possible flux and spectral variations of the INS in detail.

We used the SAS task eposcorr to refine the astrometry by cross-correlating the list of EPIC X-ray source positions with those of catalogued near-infrared (2MASS, Skrutskie et al., 2006), optical (GSC 2.3.2, Lasker et al., 2008), and X-ray (Chandra, Evans et al., 2010) objects lying within $15^{\prime}$ from J1605. The results that yielded the least astrometric errors were obtained when cross-correlating the X-ray sources with a number of around 50 optical counterparts present in the field-of-view. Small positional offsets in right ascension and declination were consistently detected for all catalogues under study and are also shown in Table 3. The astrometrically corrected EPIC source positions in all four observations are consistent with each other.

Finally, we verified the statistics of the EPIC lightcurves for general-trend variability. The lightcurves, binned into $600$ s to $1200$ s intervals, were corrected for bad pixels, dead-time, exposure, and background counts with the SAS task epiclccorr. All 2012/AO14 exposures are consistent with a constant flux.

3 Analysis and results

3.1 Timing analysis

For the timing analysis we used a $Z^{2}_{\rm m}$ test (Buccheri et al., 1983) applied directly on the times-of-arrival of the pn and EPIC (pn+MOS) events to search for periodic signals. The times-of-arrival of the pn and MOS photons were converted from the local satellite to the solar system barycentric frame using the SAS task barycen and the astrometrically corrected source coordinates in each camera and observation (Table 3).

Apart from the 2012 observation, the only other XMM-Newton datasets suitable for timing analysis are a pn LW exposure performed in 2003 (0157360401, $\sim 33$ ks) and one performed in 2002 in timing (TI) mode (0073140501, $\sim 30$ ks); these observations were used to derive an estimate of the spin down rate of the source in combination with the 2012 dataset (Pires et al., 2014). The other archival pn observations of J1605 were either operated with the thick filter (hence reducing count rates by more than a factor of two) or severely affected by background flares. Likewise, the time resolution of 2.6 s of the archival MOS observations of J1605, all performed in FF mode, is not sufficient to detect the $3.39$ s modulation.

The re-analysis of the 2012 observation with SAS 15 and up-to-date calibration files yields similar results as those reported in Pires et al. (2014). The only statistically significant modulation, at $\nu\equiv\nu_{2012}=0.2951709(14)$ Hz, results when the search is restricted to source photons with energy above 0.5 keV. We also verified that the periodic signal is always present at the same frequency within the errors, independently of the exact details of the processing of the raw event file (e.g. included calibration files, SAS version, randomisation in energy within a PI channel, event filtering, or randomisation in time within the sampling detector time). The measured fluctuation of the power of the $Z^{2}_{1}$ statistic at $\nu_{2012}$, $Z^{2}_{1}(\nu_{2012})\sim 35-49$ is consistent with that expected from a sinusoidal modulation of amplitude $p_{\rm f}=(4.1\pm 0.9)$% (see, for example, Pavlov et al., 1999).

We first analysed each of the four AO14 observations individually. Taking into account the typical spin-down rate of the M7 INSs (see, e.g. Kaplan & van Kerkwijk, 2009, and references therein), we looked for significant signals in a $5\times 10^{-4}$ Hz range around the 2012 frequency, adopting a resolution of $0.1$ $\mu$Hz (oversampling factor of at least 10). The number of statistically independent trials in each search, which depends on the frequency range and on the total duration of the observation, is typically within 30 to 60. Assuming the usual scenario of magnetic dipole braking in vacuum, the $Z^{2}_{\rm m}$ tests allow for a maximum braking with respect to the 2012 signal of $|\dot{\nu}|\lesssim 3.5\times 10^{-12}$ Hz s*-1*, which corresponds to that exerted by a maximum dipolar magnetic field of $B_{\rm dip}\lesssim 4\times 10^{14}$ G at the equator.

Due to the energy-dependent nature of the 2012 signal (see Pires et al., 2014, for details), we carried out tests in various energy bands, also varying the size of the source extraction region and other parameters of the search (e.g the included photon patterns and other details of the processing of the raw event files). We tested seven energy ranges in the soft ($0.15-0.5$ keV, $0.2-0.5$ keV, and $0.3-0.5$ keV), hard ($0.5-1.35$ keV), and total ($0.15-1.35$ keV, $0.2-1.35$ keV, and $0.3-1.35$ keV) energy bands; source counts – between $(0.18-5.5)\times 10^{5}$ pn and $(0.4-7.2)\times 10^{5}$ EPIC photons – were extracted from circular regions of radius $10^{\prime\prime}$, $25^{\prime\prime}$, $50^{\prime\prime}$, $100^{\prime\prime}$ and $120^{\prime\prime}$ around the source position in each camera and epoch. Altogether, 1680 tests were conducted in the pn/EPIC datasets of the four observations.

The summary of results is in Figure 1. In each plot we show the $Z^{2}_{1}$ statistic as a function of trial frequency in each AO14 observation (pn/EPIC datasets), for different energy bands, and taking the extraction region of radius $100^{\prime\prime}$ as an illustrative example. The $Z^{2}$-test performed on the 2012 pn observation and respective peak at $Z^{2}_{1}(\nu_{2012})\sim 50$ (out of scale) is shown in the background of the four plots for comparison ($0.5-1.35$ keV). No significant power (above $3\sigma-4\sigma$) in the searched frequency range is consistently detected in these tests. The inclusion of higher harmonics (or different photon patterns and event filters) in the $Z^{2}_{\rm m}$ tests does not affect the results. Therefore, unless significant changes of pulsed fraction have taken place, $4\sigma$ upper limits666See, e.g. Groth (1975), for the method to extract upper limits on the pulsed fraction. of $1.33(6)$%, $1.74(8)$%, $1.84(8)$%, and $1.80(8)$% in the total energy band in each AO14 observation rule out the 2012 candidate period at a high confidence level.

We next searched the new data, especially the longest and most sensitive ‘201’ observation, for other significant modulations in the full frequency range allowed by the resolution of the EPIC cameras (‘blind searches’). To this end, the times-of-arrival of the pn and MOS events ($\sim 2.5\times 10^{5}$ counts) were analysed together in the $\nu=0.01-0.56$ Hz frequency range; for pn, timing searches were extended to higher frequencies, up to $\sim 6.8$ Hz888The maximum frequency is determined by the Nyquist limit, given the time resolution of the detector.. The adopted frequency step was $\Delta\nu=2.5-5$ $\mu$Hz (oversampling factor of 3) and the number of independent trials were $(4-8)\times 10^{5}$ and $(3-6)\times 10^{4}$ in the pn and EPIC searches, respectively. A total of seven ‘narrow’ (with widths between $100$ eV and $600$ eV) and three ‘wide’ ($600-850$ eV) energy bands, defined within the $0.15-1.35$ keV range according to the source’s signal-to-noise ratio, were defined for these searches. We tested three different sizes of source extraction regions ($25^{\prime\prime}$, $50^{\prime\prime}$, and $120^{\prime\prime}$) and included all valid photon patterns in the EPIC searches. For pn, we restricted the event lists to include only single and double photon patterns. In total, 240 EPIC and pn blind searches were performed in the observations of the large programme.

No significant periodic signal resulted from the analysis. In Table 4 we list the most constraining $3\sigma$ upper limits from the EPIC and pn searches for three wide energy band intervals (soft, hard, and total; see the table caption for details).

3.2 Spectral analysis

3.2.1 EPIC data

The analysis of the EPIC data is based on source and background spectra extracted from regions as described in Section 2, together with the respective response matrices and ancillary files created for each of the EPIC cameras and observation. In accordance to the guidelines and calibration status of the instruments, we restricted the spectral analysis to GTI-filtered photons with energies between 0.3 keV and 1.35 keV (beyond which the source signal-to-noise ratio becomes insignificant). Including the 2012 observation, the analysed dataset comprises 15 spectra and over $1.2\times 10^{6}$ counts ($0.3-1.35$ keV), of which 1.8% can be ascribed to the background.

The pile-up level is negligible in the MOS exposures (Section 2); for pn, we applied a correction in the redistribution matrix files with the SAS task rmfgen to minimise flux loss and spectral distortion999https://www.cosmos.esa.int/web/xmm-newton/

sas-thread-epatplot. The correction is calculated directly from the frequency and spectrum of the incoming photons and has the advantage of keeping events from the central PSF area in the spectral analysis, which would have to be otherwise discarded. We verified that the results using this approach agree well with those when the spectra are extracted from regions with a $10^{\prime\prime}-15^{\prime\prime}$ excised core, while avoiding up to $30\%$ of data loss in the pn exposures.

The energy channels of each spectrum, which are by construction 5 eV wide, were regrouped to avoid a low ($<30$) number of counts per spectral bin. Due to the brightness of the source and the good statistics of each individual spectrum this has an effect only at the high energy side of the analysis, where the source signal becomes dominated by the background. At lower energies ($\lesssim 700$ eV), the spectrum oversamples the instrument resolution of the EPIC cameras by a factor of up to 20. We ensured nonetheless that oversampling did not influence the results of spectral fitting – specifically, we checked for consistency where oversampling was kept within a maximum factor of 3.

To fit the spectra we used XSPEC 12.9.0n (Arnaud, 1996). Unless otherwise noted, the fit parameters were allowed to vary freely within reasonable ranges. The photoelectric absorption model and elemental abundances of Wilms et al. (2000, tbabs in XSPEC) were adopted to account for the interstellar material in the line-of-sight. Due to the low absorption towards J1605, the choice of abundance table and cross-section model does not significantly impact the results of the spectral fitting.

The 15 spectra were first fit individually to check the agreement between the instruments and epochs. The exercise showed the expected few percent cross-calibration uncertainty between the EPIC detectors (Read et al., 2014). Although variations from pointing to pointing for a given instrument are formally significant with respect to a constant value, the relative error is still smaller than the absolute discrepancies between the cameras within a given epoch. To account for this uncertainty, we allowed for a renormalisation factor in XSPEC and fitted the spectra of the pn and MOS cameras simultaneously. We checked that the results of the simultaneous fits were consistent with the weighted means of the individual measurements.

We then checked if consistent results are obtained when fitting a single ‘stacked’ spectrum, which converges to the best parameter values much faster in XSPEC than the simultaneous fits. The stacked spectra are produced with the SAS task epicspeccombine, taking into account the responses, background, and ancillary files of the individual exposures. For pn, the stacking approach leads to inconsistent results that do not match the corresponding weighted mean values for the camera within the errors, nor the results of the simultaneous fits111111See https://www.cosmos.esa.int/web/xmm-newton/

sas-thread-epic-merging#cav, for details.. On the other hand, the results of the combined MOS spectra agree well with those of the simultaneous fits. We adopted hereafter a stacked spectrum only for MOS to avoid introducing biased results in the spectral analysis.

Next, we proceeded at finding a model which closely describes the X-ray spectral energy distribution of the source (see Pires et al. 2014, for details). We list in Table 5 the results of the fit of a double blackbody model131313For comparison, we also list in Table 5 for the multi-epoch fits the results of a single temperature blackbody model with a Gaussian absorption line. For this model, the fit quality is generally poor ($\chi^{2}_{\nu}\gtrsim 2$) and the column density is unconstrained. with a broad Gaussian absorption feature. In the fitting procedure we restricted the energy of the line between $0.3$ keV and $1.35$ keV and its Gaussian $\sigma$ between [math] eV and $200$ eV; the column density is varied between $N_{\rm H}=0$ cm*-2* and $5\times 10^{21}$ cm*-2*, while the temperature of the blackbody components can assume values between $5$ eV and $500$ eV. For each observation, we fitted the model to the pn spectrum, to the combined MOS1 and MOS2 spectrum (labeled ‘MOS’ in Table 5 and the two subsequent tables), and to the pn and MOS spectra simultaneously (labeled ‘EPIC’ in Table 5 and the two subsequent tables). Finally, we performed multi-epoch simultaneous fits of all pn (5 spectra), MOS1/2 (10 spectra), and EPIC (6 spectra, comprising 5 pn and one combined MOS) exposures.

For each fit in Table 5, we list the reduced chi-square ($\chi^{2}_{\nu}$) and its null-hypothesis probability (NHP in %), the column density $N_{\rm H}$ in units of $10^{20}$ cm*-2*, the temperature of the cold $kT_{1}^{\infty}$ and hot $kT_{2}^{\infty}$ blackbody components in eV, the radiation radii $R_{1}^{\infty}$ and $R_{2}^{\infty}$ of each component (assuming a distance to the source of $d\equiv d_{300}=300$ pc; Tetzlaff et al. 2012, e.g.), the central energy of the absorption line (when constrained) or the corresponding $3\sigma$ upper limits, its Gaussian $\sigma$, and equivalent width $EW$ (all in eV), and the observed model flux in the energy band $0.2-12$ keV, $f_{\rm X}$, in units of $10^{-12}$ erg s*-1* cm*-2*. The model provides acceptable $\chi^{2}_{\nu}$ values in each epoch, with NHP between 14% and 88%. The somewhat large chi-square values and worse fit quality (NHP $<1\%$) of the multi-epoch fits could not be improved by the inclusion of an additional model component.

In the left column of Figure 2 we plot the results of Table 5 as a function of time, with $1\sigma$ errors. The best-fit parameters per instrument are consistent between pointings despite the systematic differences between the detectors. The column density is constrained and for the MOS and EPIC fits agrees within errors with the Galactic value in the direction of the source, $N_{\rm H}^{\rm gal}=(2.4-2.6)\times 10^{20}$ cm*-2* (e.g. Kalberla et al., 2005; Willingale et al., 2013). The pn camera measures twice as much absorption, as well as $7\%$ to $9\%$ softer temperatures with respect to MOS; the observed model flux is also $6\%$ lower. If $N_{\rm H}$ is fixed to the Galactic value, the temperature of the two components and the observed model flux typically agree within $2.5\%$. Considering the simultaneous EPIC fit as an effective average between the detectors, the best-fit parameters are well constrained within ranges $N_{\rm H}=(2.4-4)\times 10^{20}$ cm*-2*, $kT_{1}^{\infty}=(64-68)$ eV, $kT_{2}^{\infty}=(120-124)$ eV, $R_{1}^{\infty}=(9.5-13)$ km, $R_{2}^{\infty}=(1.0-1.2)$ km, $\epsilon<(340-420)$ eV, $\sigma=(100-120)$ eV, $EW<(150-200)$ eV, and $f_{X}=(6.4-6.9)\times 10^{-12}$ erg s*-1* cm*-2*. The model flux corrected for absorption (unabsorbed) is $F_{\rm X}=(1.2-1.7)\times 10^{-11}$ erg s*-1* cm*-2*.

Alternatively, the dataset can be as well fit (NHP $\sim 10\%-80\%$) by a fully ionized neutron star hydrogen atmosphere model (nsa in XSPEC; Pavlov et al., 1995; Zavlin et al., 1996), again modified by a broad Gaussian absorption line (Table 6). We tested non-magnetised ($B<10^{8}$ G) and magnetised models with magnetic field values of $B=10^{12}$ G and $B=10^{13}$ G; in the fitting procedure, the neutron star mass and radius were at first fixed at the canonical values, $M=1.4$ M⊙ and $R=10$ km, and then allowed to vary to check for an improved fit. In Table 6 we list the results for a canonical neutron star with $B=10^{13}$ G, which is the model that in most cases gave the highest NHP for each dataset. The unredshifted model effective temperature $T_{\rm eff}$ in K, the distance $d$ in pc, the parameters of the line ($\epsilon$, $\sigma$, $EW$), and the observed model flux $f_{\rm X}$, are also listed in Table 6.

In the right column of Figure 2 we plot the best-fit nsa parameters as a function of MJD (assuming for all epochs the results of the $B=10^{13}\,G$ fits with $M=1.4$ M⊙ and $R=10$ km). In comparison with the double blackbody model, the systematic differences between the detectors persist, however the measurements differ by a much smaller percentage – $13\%$ in $N_{\rm H}$, $\lesssim 2\%$ in $T_{\rm eff}$, and $3\%$ in $f_{\rm X}$ – and the overall consistency between epochs and instruments is improved (c.f. overlapping shaded areas). The central energy of the line is well constrained, and also narrower than in the double blackbody case. Again considering the best-fit results of the EPIC fits we have: $N_{\rm H}=(2-3)\times 10^{20}$ cm*-2*, $T_{\rm eff}=(5.7-5.8)\times 10^{5}$ K, $B=10^{13}$ G, $d=(124-129)$ pc, $\epsilon=(360-400)$ eV, $\sigma=(85-100)$ eV, $EW=(55-60)$ eV, and $f_{X}=(6.5-7.0)\times 10^{-12}$ erg s*-1* cm*-2*. The unabsorbed flux of this model is measured in a similar range as that of the double blackbody model, $F_{\rm X}=(1.2-1.4)\times 10^{-11}$ erg s*-1* cm*-2*. In Figure 3 we show this best-fit model folded to the EPIC dataset with residuals.

The neutron star distance derived from the nsa fits, $d\sim 110-130$ pc, is rather small in comparison to the range expected for the source, $d\sim 300-400$ pc (Posselt et al., 2007; Tetzlaff et al., 2012). While the fit is insensitive to the mass of the neutron star, unrealistically large neutron star radii, $R>20$ km, and an overall poor fit quality, $\chi^{2}_{\nu}\sim 1.7$, are obtained when the distance to the source is fixed at around the estimated value (see discussion in Section 4).

We explored other neutron star models where the atmosphere can be partially ionized and the size of the emission radius, $R_{\rm em}$, can be parametrised with respect to the neutron star physical radius (nsmaxg in XSPEC; Ho et al., 2008). For each epoch and instrument, as well as for the multi-epoch fits as before, we tested 19 absorbed nsmaxg models with $B=(0.01-30)\times 10^{12}$ G and $M=1.4$ M⊙. We first set the size of the emitting region to be the same as the neutron star radius; then we allowed this parameter to vary to smaller values to check for improved fits. As for the nsa models, we show in Table 7 the results with $B=10^{13}$ G, which are the ones with the generally highest NHP and the most consistent parameters between epochs and instruments.

The best-fit models (with somewhat comparable NHP as in the 2bb and nsa models, i.e., between $2\%$ and $89$%) are for a neutron star atmosphere composed of hydrogen at a distance of less than $110$ pc ($3\sigma$). All models consisting of mid-$Z$ element plasma (C, O, Ne) provided poor fit results. While the properties of the absorption line were found to be nearly identical to those of the fully ionized case, the temperature of the atmosphere is higher, and the radiation is roughly twice as much absorbed and inconsistent with the Galactic value. The size of the emission region was found to be slightly smaller than the canonical 10 km of the nsa models, with $R_{\rm em}\sim 8-9$ km. The best-fit results of the EPIC fits are within: $N_{\rm H}=(3-4)\times 10^{20}$ cm*-2*, $T_{\rm eff}=(6.0-6.3)\times 10^{5}$ K, $B=10^{13}$ G, $d<110$ pc, $\epsilon=(340-400)$ eV, $\sigma=(85-100)$ eV, $EW=(50-65)$ eV, and $f_{X}=(6.4-6.8)\times 10^{-12}$ erg s*-1* cm*-2*. The unabsorbed flux, $F_{\rm X}=(1.3-1.5)\times 10^{-11}$ erg s*-1* cm*-2*, is consistent with those of the other two previously discussed models.

To break some of the degeneracy between the parameters and look for more physical results, we restricted the distance to the source within $d=100-600$ pc, capped the column density at the Galactic value, and let the neutron star mass and radius vary within $M=0.5-2.5$ M*⊙* and $R=5-15$ km. However, the exercise led to generally worse fits and significant discrepancies between the best-fit parameters of pn and MOS. No other neutron star atmosphere model in XSPEC, nor the inclusion of a second (colder) component, provided acceptable fits.

3.2.2 RGS data

For the spectral analysis of RGS data we included four XMM-Newton observations of the source performed in 2002/2003 in addition to the five 2012/AO14 observations, thus considerably extending the time span of the analysis in relation to that covered by the EPIC data (Tables 1 and 2 in Section 2). The total analysed RGS dataset, of which the AO14 observations account for nearly 70% in net exposure, amount to 18 RGS1/2 spectra and GTI-filtered exposures of 465 ks and 457 ks per detector.

We used the EPIC source coordinates in each observation to generate the instrument spatial masks and energy filters with rgsproc. The GTI-filtered event lists were used to extract the source and background spectra in wavelength space using the tasks rgsregions and rgsspectrum, while response matrix files were produced with the SAS task rmfgen. Only the first-order spectra were analysed. To increase the signal-to-noise ratio each spectrum was rebinned into 0.165 Å wavelength channels. The defective channels of the RGS cameras141414XMM-Newton Calibration Technical Note 0030, issue 7.7; hereafter Gonzàlez-Riestra et al. (2018)., which cover in first order the wavelength ranges of $11$ Å to $14$ Å in RGS1 and $20$ Å to $24$ Å in RGS2, were excluded from the spectral fitting. The total dataset amount to $3.992(20)\times 10^{4}$ and $3.724(19)\times 10^{4}$ counts ($15-30$ Å), respectively in each RGS1/2 camera, of which around $40\%$ can be ascribed to the background.

We fitted each observation in XSPEC assuming a model (hereafter, the bbgauss model) consisting of an absorbed blackbody, modified by a Gaussian absorption feature with $\sigma=100$ eV as found from the analysis of EPIC data; the column density was fixed to the Galactic value to better constrain the other model parameters. The RGS1/2 spectra were fitted simultaneously adopting a constant factor between the instruments. We note that the best-fit parameters from the fit of individual RGS1/2 spectrum agree well with each other in a given epoch.

The best-fit results of the bbgauss model are in Table 8. For each observation, we list the blackbody temperature $kT_{\infty}$ in eV, the radiation radius $R_{\infty}$ in km (assuming a source at $d_{300}$), the central energy $\epsilon$ and equivalent width $EW$ of the Gaussian absorption feature, and the unabsorbed source flux of the model $F_{\rm X}$, in the $0.2-12$ keV energy band. The results suggest a possible trend of the parameters of the source: in comparison with the first four observations obtained in 2002 and 2003 and labeled (A-D) in Table 8, the more recent observations show a slight increase in temperature and a more pronounced decrease in the model normalization and flux (formally inconsistent with a constant term), at constant properties of the Gaussian absorption feature (Figure 4). By contrast, the parameters of the source within these two subgroups are constant at the $2\%-3\%$ level in $kT_{\infty}$, $7\%-10\%$ in $R_{\infty}$, and $5\%-10\%$ in $F_{\rm X}$. The trends are seen in the spectral parameters of both RGS1/2 instruments.

In Pires et al. (2014), we investigated the constancy of the INS emission on the EPIC data performed between 2002 and 2012. Unfortunately, the analysis does not allow us to draw definite conclusions: while the MOS instruments are unsuited for long-term studies161616XMM-Newton Calibration Technical Note 0018., the pn camera provides only one data point prior to 2012 for comparison, due to the heterogeneous dataset and background flares (Sections 1 and 2). Nonetheless, an increase of blackbody temperature, consistent with what is observed in the RGS data, was then reported.

To investigate the possibility of a long-term evolution on the parameters of J1605, we co-added the spectra of the two subgroups in each RGS camera, using the SAS task rgscombine, and binned the results to 0.165 Å as before. The resulting grouped spectra (labeled ‘old’ and ‘new’ in Table 8) were then fitted simultaneously in XSPEC with the same bbgauss model, taking into account the co-added background and response files in each detector as usual. In Figure 5 we plot the ‘old’ and ‘new’ spectra (grey and magenta data points) of J1605, with the folded bbgauss model and fit residuals. The best-fit results as a function of time are plotted as shaded grey and pink areas in Figure 4.

The results of this approach confirm the observed trend. With respect to the early pointings, we measure a $7\%$ higher temperature and a $25\%$ lower flux and smaller radiation radius in the observations performed after 2012, formally significant beyond the spectral errors. Nonetheless the significance of the variations is low: $1\sigma$ in $kT$ and $F_{X}$, while the other parameters are consistent within the rather large spectral errors (e.g. the emission radius in Figure 4).

The RGS instruments suffer from a decline in sensitivity at long wavelengths, likely due to a build-up of of hydrocarbon contamination on the detector (de Vries et al. 2015, Gonzàlez-Riestra et al. 2018). Empirical corrections were first introduced in 2006 to take this and other effects into account in the calibrated model of the RGS1/2 effective areas, which are estimated to have an absolute accuracy of 10%. Indeed, we observe a 18% decrease in sensitivity in the effective area at long wavelengths between the old and new datasets (apparent in Figure 5). Altogether, uncertainties from both the spectral model and other calibration issues, possibly not accounted for in the modeling of the effective area with time, may be responsible for the discrepancies on the parameters of J1605 reported here.

The inclusion of a cold blackbody component in the bbgauss model, unlike for the EPIC data, is not satisfactory due to the large normalization required to fit the RGS spectra. If this is kept within reasonable limits (that is, corresponding to a $<10^{33}$ erg s*-1* blackbody at $d=0.1-1$ kpc), the quality of the fit is worsened; moreover, while there are no significant changes on the temperature of the hot component and on the parameters of the absorption line, the best-fit temperature of the cold component is very soft, $<30$ eV, and the column density is 2 to 3 times higher than the Galactic value.

The evidence for a narrow absorption feature at energy $\epsilon\sim 0.57$ keV ($\lambda=21.5$ Å) in the RGS spectra of J1605 was first reported by van Kerkwijk et al. (2004). Similarly narrow features at around this wavelength have been identified in the RGS spectrum of the M7 INS RX~J0720.4-3125 (Hambaryan et al., 2009) and other thermally emitting INSs (Hohle et al., 2012). To investigate the presence of the narrow feature in our dataset, we fitted each of the nine RGS1171717The RGS2 data cannot be used due to the defective channels of the camera around the wavelength range of interest. spectra individually, using an absorbed blackbody model and two Gaussian absorption lines. As the individual datasets do not have very high signal-to-noise, we fixed the column density to the Galactic value and set the energy and FWHM of the broad absorption line to the best parameters found consistently in the analysis of EPIC and RGS data ($\epsilon=410$ eV and $\sigma=100$ eV). Absorption features were then searched between 21 Å and 22.5 Å ($550-590$ eV).

The results are summarised in Table 9. For each fit, we show the signal-to-noise $S/N$ of each spectrum, the blackbody temperature $kT$ of the source and the equivalent width of the broad absorption feature $EW_{1}$; the energy $\epsilon_{2}$, FWHM $\sigma_{2}$, and equivalent width $EW_{2}$ of the narrow feature – when constrained, or their corresponding $1\sigma$ upper limits – are also listed.

The best-fit spectral parameters are consistent with those of Table 8, showing that the inclusion of the narrow feature is not statistically required in most cases. Considering the number of trials ($40$) in the searched wavelength range, the evidence for the narrow feature is only significant in observation (A). Remarkably, there is no evidence for a narrow feature within $21-22.5$ Å in the longest observation of the source (labeled ‘201’), which has a much higher $S/N$ than the others. The same analysis carried out in the co-added ‘old’ and ‘new’ spectra confirm that additional features are absent in the most recent pointings.

4 Discussion

The M7 have been considered a rather homogeneous group of cooling neutron stars, displaying similar ages, temperatures, and timing properties. The source RX~J1605.3+3249 stood out in that it could be slowing down at a fast rate, indicating a high dipolar field – the highest amongst the group – and a possible evolution from a magnetar (Pires et al., 2014). The analysis of our dedicated XMM-Newton large programme does not confirm the previous results (Section 3.1). Due to the energy-dependent nature of the previously detected modulation, we performed extensive high-resolution periodicity searches allowing for moderate changes of pulsed fraction and the optimal energy range and signal-to-noise ratio for detection, for a reasonably wide range of spin down values. No significant signal resulted from the analysis: unless considerable changes of pulsed fraction have taken place since 2012, the deepest upper limit of $1.33(6)\%$ ($4\sigma$), in the relevant frequency range, conservatively rules out the $3.39$ s modulation. Moreover, in the full frequency range allowed by the timing resolution of the EPIC cameras, blind searches revealed no other periodic signals with $p_{\rm f}\gtrsim 1.5\%$ ($3\sigma$; $0.3-1.35$ keV), thus considerably improving previous estimates for pulsations with $P>0.15$ s. Similarly low $3\sigma$ upper limits, within $1.8\%$ and $4\%$, are obtained in the same period range in narrow (100 eV to 600 eV wide) energy intervals, defined according to the source’s signal-to-noise ratio.

With over $10^{6}$ EPIC counts, the unprecedented photon statistics of the new dataset allowed the deepest to-date investigation of the X-ray emission of the source (Section 3.2.1). We found that, altogether, no theoretical model available in XSPEC can provide a fully satisfactory physical description of the spectrum of the neutron star. While statistically acceptable fits are obtained for the individual epochs, multi-epoch fits including the spectra of all EPIC cameras have null-hypothesis probabilities of less than $1\%$, which may at least partially be ascribed to cross-calibration uncertainties. Best results were obtained when fitting the data with either a double-blackbody (2bb) or a magnetised neutron star atmosphere model consisting of hydrogen (nsa and nsmaxg, with $B=10^{13}$ G and $M=1.4$ M⊙), in either case modified by a broad Gaussian absorption feature as previously reported in the literature.

No significant evidence of spectral variability is measured in the 2012–2016 time frame covered by the analysis of EPIC data. The overall consistency of the parameters (between epochs and EPIC instruments) was optimal for the atmosphere models, with systematic errors of $13\%$ in column density, $2\%$ in temperature, and $3\%$ in flux. In particular, a canonical nsa model with $B=10^{13}$ G much better constrains the column density toward the source and the properties of the absorption line in relation to the 2bb model. Best-fit distances around $d\sim 130$ pc are, on the other hand, inconsistently smaller than that estimated for the source, $300\pm 50$ pc. The model consisting of a partially ionized hydrogen neutron star atmosphere (nsmaxg) provides even smaller distances ($3\sigma$ upper limits below $110$ pc) for a $R_{\rm em}\sim(8-9)$ km emitting region on the neutron star, while the derived column density is nearly two times the Galactic value in the direction of the source.

Considering the nsa model of the multi-epoch pn fits (which is the result with the highest NHP among the multi-epoch fits), the spectral parameters of the source are constrained as $N_{\rm H}=2.45(14)\times 10^{20}$ cm*-2*, $T_{\rm eff}=5.81(3)\times 10^{5}$ K, $\epsilon=385\pm 10$ eV, $\sigma=86.8\pm 0.3$ eV, $EW=55\pm 3$ eV, and $f_{\rm X}=6.842(9)\times 10^{-12}$ erg s*-1* cm*-2* ($0.2-12$ keV). In contrast to previous analysis, we did not find that the inclusion of other model components (in particular, additional lines in absorption) were statistically justified or could significantly improve the results of the multi-epoch fits. Other up-to-date, fully and partially ionized neutron star atmosphere models, consisting of different elemental compositions and with non-canonical neutron star mass and radius, did not provide better fits than the models described above.

The typical distance derived from the best-fit atmosphere models, around 130 pc, is smaller than the range expected for the source, $d=350\pm 50$ pc (Posselt et al., 2007, Section 3.2.1). This range is derived from the fitted hydrogen column density, assuming a blackbody model with three Gaussian lines in absorption, and a three-dimensional description of the distribution of the interstellar medium in the direction of the source (which, according to the authors, should be reliable up to $\sim 270$ pc). Based on kinematic arguments, Tetzlaff et al. (2012) applied this expected range and the observed proper motion of the source (Motch et al., 2005; Zane et al., 2006) to trace back the neutron star trajectory and determine its likely birthplace, using possible associations with runaway massive stars and the observed abundance of heavy elements as further evidence to narrow down the most likely solutions. These predict a current distance of $d=300-370$ pc if the neutron star was born less than 0.5 Myr ago in a nearby supernova explosion.

4.1 The viewing geometry and presence of hotspots

Magnetic fields in the range observed in the M7 are expected to produce large temperature variations on the neutron star surface, due to the anisotropic electron conductivity and heat transport in the stellar envelope and crust (e.g. Geppert et al., 2004; Pérez-Azorín et al., 2006; Perna et al., 2013). In this case, strong pulsed flux variations are expected at the neutron star spin period unless the source is observed from a particularly unfavourable geometry: either if the angle between the neutron star spin axis and the line-of-sight $i$ is sufficiently small, or if the regions of higher temperature (hotspots) are located at a very small angle $\theta_{B}$ in relation to the neutron star rotation axis (that is, the magnetic and spin axes of the star are nearly co-aligned). Therefore, the stringent pulsed fraction limits from the timing analysis can be used to verify the viability of the 2bb model and probe the presence of hotspots on the surface of J1605 (e.g., see the case of the thermally emitting central neutron star in the supernova remnant HESS~J1731-347, Suleimanov et al., 2017).

With this goal we considered the model of a slowly rotating neutron star, observed at an inclination angle $i$, with two identical polar hotspots located at an angle $\theta_{B}$ with respect to the spin axis (see, e.g. Page, 1995; Schwope et al., 2005; Suleimanov et al., 2010, for a full description of the model). Light bending in the vicinity of the neutron star follows the relation between the local angle of the emitted photon and its escape direction and depends on the compactness of the neutron star (given by the ratio between the neutron star and the Schwarzschild radius $r_{\rm g}=R_{\rm ns}c^{2}(2GM_{\rm ns})^{-1}$; Beloborodov 2002). The temperature of the neutron star surface, $kT_{\rm ns}^{\infty}=60.9_{-1.5}^{+1.7}$ eV, and of the hotspots, $kT_{\rm sp}^{\infty}=117.0\pm 0.8$ eV, are assumed from the best-fit 2bb pn model, which is the one with the highest NHP (Table 5). For a neutron star distance of $d=300^{+50}_{-40}$ pc (Tetzlaff et al., 2012), the model normalisations set the corresponding sizes of the emission regions as $R_{\rm ns}^{\infty}=16.4_{-2.2}^{+2.9}$ km and $R_{\rm sp}^{\infty}=1.35_{-0.18}^{+0.24}$ km. With these values the angular size of the spots and the compactness of the star are fixed as $\theta_{\rm sp}=4.7^{\circ}$ and $r_{g}=2.8$, respectively, assuming a 1.5 M⊙ neutron star.

For a particular viewing geometry ($i$, $\theta_{B}$), the photon flux at a given rotation phase results from the sum of the visible individual area elements of the neutron star surface, assuming blackbody emission and taking into account the light bending. The flux is corrected for the interstellar absorption of a equivalent column density of $N_{\rm H}^{\sf 2bb}=4.5\times 10^{20}$ cm*-2* and then folded with the EPIC pn response to give the source count rate at the $0.5-1.35$ keV energy band. The energy band is chosen to minimise the effects of the broad absorption feature in the emitted spectrum.

With this method we computed an extensive grid of lightcurves for $(\theta_{B},i)$ within $(0^{\circ},0^{\circ})$ and $(90^{\circ},90^{\circ})$ and computed the pulsed fraction for each orientation as:

[TABLE]

In Figure 6 we plot the resulting $p_{\rm f}$ map in the $(\theta_{B},i)$ plane. The maximum pulsed fraction obtained for the 2bb model is about 20%. The region allowed by the limits of the timing analysis lies below the thick pink line corresponding to $p_{\rm f}=2.78(16)$% (Table 4; $0.5-1.35$ keV). Integrating over all possible random orientations of line-of-sight inclination and spot angles we obtain a small likelihood ($\sim$1.9%) that we do not see pulsations from the source due to the particularly unfavourable viewing geometry, if the 2bb model is correct.

4.2 The energy distribution

Thermal emission from INSs is expected to originate immediately at the surface, with the bulk of the energy flux peaking in the soft X-ray band. In principle, by confronting the observed spectra and light curves with theoretical models for neutron star thermal radiation, it should be possible to derive the surface temperature, magnetic field, gravitational acceleration and chemical composition: if distances are known, then the stellar mass, radius and the equation of state of neutron star interior could be constrained as well (see Potekhin et al., 2015; Özel & Freire, 2016, for recent reviews on neutron star atmosphere models and up-to-date astrophysical constraints on the equation of state of nuclear matter). Since their discovery in the All-Sky Survey of the ROSAT satellite (Voges et al., 1999), the M7 have been regarded as the closest-to-perfect candidates for testing neutron star emission models, due to a combination of bright thermal emission, proximity, independent distance estimates191919From HST parallaxes (in two cases) and kinematic studies (e.g. Walter et al., 2010; Tetzlaff et al., 2010, and references therein)., and a lack of significant magnetospheric or accretion activity.

In practice, progress has been hampered by uncertainties on the chemical composition of the atmosphere and the phase state of the stellar surface, as well as by the lack of understanding on the magnetic field and temperature distributions (e.g. Zane & Turolla, 2006; van Kerkwijk & Kaplan, 2007; Suleimanov et al., 2010). The presence of lines adds to this complexity as, although believed to be related with the star’s magnetic field, they have no unique physical interpretation. To explain the emitted radiation and equivalent widths of the lines in the phase-resolved spectrum of the M7 RXJ1308.6+2127, Suleimanov et al. (2010) favoured a model where a partially ionized, optically thin atmosphere above the condensed surface must be present (see also Motch et al., 2003; Ho et al., 2007). Using this model, Hambaryan et al. (2011) derived the temperatures of the X-ray emitting areas and the magnetic field intensity at the poles; moreover, they could constrain the compactness of the neutron star and the gravitational redshift on the surface, suggesting a very stiff equation of state. Similar conclusions were reached for the M7 RXJ0720.4-3125 (Hambaryan et al., 2017). These results are only marginally compatible with the most favoured range of the true radius of a 1.5 M⊙ neutron star, $10-11.5$ km, from the analysis of Özel & Freire (2016). New generation X-ray missions, in particular the Neutron Star Interior Composition Explorer Mission (NICER, Gendreau et al., 2012), together with more accurate distances from the GAIA satellite (Gaia Collaboration et al., 2016), will certainly improve the constraints on neutron star mass and radius from astrophysical observations of, for example, millisecond pulsars in globular clusters.

Recently, Viganò et al. (2014) showed that in some cases (as for the M7 RX~J0806.4-4123; Haberl et al. 2004) the deviations found in the spectra of thermally emitting INSs may be induced simply by the inhomogeneous temperature distribution on the surface. While the effect is unlikely to account for all cases of sources with reported spectral features, the interesting result is that the anisotropic temperature distribution can give way to “spurious” spectral features to be claimed. We can safely exclude this possibility for the absorption feature in J1605, which cannot be accommodated by a multi-temperature energy distribution.

All M7 INSs have detected optical, ultraviolet, or infrared counterparts (see Kaplan et al., 2011; Posselt et al., 2014, 2018, for references and limits). Interestingly, the extrapolation to longer wavelengths of the best-fit model inferred from X-rays – including, for the case of J1605, both the double temperature blackbody and atmosphere models discussed here – falls below the actual detected fluxes: this is known as the ‘optical excess’, and is observed in all M7 INSs. The optical excess of J1605 deviates significantly from the expected Rayleigh-Jeans slope of the spectra and can be described by a rather flat power-law (Kaplan et al., 2011). The origin of the excess flux might rely on atmospheric effects, magnetospheric emission, or resonant scattering. In particular, if the X-ray and optical/UV emission came from different regions on the surface (e.g. Braje & Romani, 2002), we might expect correlations between the amount of optical excess and the X-ray pulsed fraction, which are not verified. The possible presence of fossil fallback disks, of a faint pulsar wind nebula, ‘bare’ neutron star surfaces, and other alternative scenarios remain open (see Turolla et al., 2004; Ertan et al., 2017; Wang et al., 2017; Posselt et al., 2018, for a discussion).

4.3 Narrow absorption feature

In combination with the existing archival RGS data, the AO14 campaign accumulated over $450$ ks of exposure time on J1605, increasing the available data by $70\%$. The good statistics allowed a detailed analysis of the narrow feature at $\epsilon\sim 570$ eV, previously reported in the literature. The investigation shows that the feature is only significantly detected in one early epoch (the 2002 observation first analysed by van Kerkwijk et al. 2004), or when this observation is combined with the archival data obtained prior to the AO14 campaign (Hohle et al., 2012; Pires et al., 2014). The feature is definitely not present in the AO14 observations, while evidence of a less significant narrow feature is present at energy $\epsilon=550$ eV in the 2012 observation. Consistently with these results, the analysis of the two grouped spectra of early and recent observations constrain the presence of the narrow feature only in the first subgroup (Section 3.2.2).

Hambaryan et al. (2009) discusses the possible physical interpretation of a similar feature detected in the co-added RGS spectrum of RX~J0720.4-3125, which was later confirmed by the analysis of Chandra LETG data (Hohle et al., 2012). Their analysis favours a blend of highly ionized oxygen originating in the ambient medium of the INS, possibly a high density nearby cloud which could contribute to the source’s optical excess. Nonetheless, a interstellar or atmosphere origins cannot be ruled out. At least for J1605, the narrow and transient nature of the feature disfavours an atmospheric origin.

Phase-dependent narrow absorption features have been reported in XMM-Newton observations of the M7 INSs RXJ0720.4-3125 and RXJ1308.6+2127 (Borghese et al., 2015, 2017), a work motivated by the detection of variable cyclotron lines detected in the spectra of two ‘low magnetic field’ magnetars (Tiengo et al., 2013; Rodríguez Castillo et al., 2016). These results give support for the presence of strong, confined magnetic field components close to the stellar surface and a complex field topology. In contrast, the features in the spectra of the two M7 INSs are intrinsically different in that they do not vary in energy, are detected at much lower energy, and are also seem to be stable and lasting over long timescales.

5 Summary and conclusions

We report here the results of a XMM-Newton large programme on the thermally emitting isolated neutron star RX~J1605.3+3249. The project aimed to gain a deeper understanding of the timing and spectral properties of the source, through a detailed analysis of its X-ray emission. The neutron star is of particular scientific interest as a source that could potentially bridge the evolutionary gap between the groups of nearby thermally emitting sources dubbed the ‘magnificent seven’ and the young and energetic magnetars. Due to the lack of detected pulsations, our science goals could only be partially completed. Nonetheless, the deep upper limits derived from our analysis were used to put stringent constraints on the viewing geometry and the presence of hot spots on the surface. Detailed phase-averaged medium and high-resolution spectroscopy constrains atmosphere neutron star models and the properties of the cyclotron line in the spectrum of the neutron star with unprecedented statistics. The non-detection of the narrow absorption feature at $\epsilon=570$ eV reported in previous epochs also disfavours an atmospheric origin.

Acknowledgements.

We thank the anonymous referee for useful comments and suggestions which helped to improve the paper. We thank Norbert Schartel and the staff of the XMM-Newton Science Operation Center, in particular Jan-Uwe Ness, for their great support in the scheduling of these time-constrained observations. We acknowledge the use of the ATNF Pulsar Catalogue (http://www.atnf.csiro.au/people/pulsar/psrcat). This work was supported by the Deutsches Zentrum für Luft- und Raumfahrt (DLR) under grant 50 OR 1511.