THE TALE OF THE TWO TAILS OF THE OLDISH PSR J2055+2539

The Large Area Telescope on board the Fermi satellite (Fermi-LAT; Atwood et al. 2009) has opened a new era for pulsar astronomy by detecting γ-ray pulsations (at E > 100 MeV) from more than 160 pulsars,⁵ ∼25% of which are not detected at radio wavelengths. The discovery of such a large sample of high-energy-emitting pulsars has enabled numerous population studies (see, e.g., Watters & Romani 2011; Pierbattista et al. 2012). This has resulted in a better understanding of the physical phenomena behind the multiwavelength emission from pulsars.

Some of the Fermi pulsars deserved dedicated multiwavelength campaigns. Such "extreme" objects, accounting for the tails of the population distribution in energetics, age, and magnetic field, are key to our understanding of the entire population, posing important limits on the different pulsar emission models. These campaigns have sometimes resulted in the discovery of important and unexpected features: e.g., the peculiar multiwavelength light-curve profile of the energetic radio-quiet pulsar PSR J1813–1246 disagrees with the classical model for X-ray emission of pulsars (Marelli et al. 2014b; Kuiper & Hermsen 2015); the isolated PSR J2021+4026 shows hints of long-term variations in the γ-ray flux, challenging current high-energy emission models (Allafort et al. 2013).

The increase in the number of γ-ray pulsars has been crucial to the study of phenomena related to highly energetic pulsars, such as X-ray and radio nebulae. Such features are usually interpreted within the framework of bow-shock, ram-pressure-dominated, pulsar wind nebulae (PWNe; for a review see Gaensler & Slane 2006). If the pulsar moves supersonically, shocked pulsar wind is expected to flow in an elongated region downstream of the termination shock (basically, the cavity in the interstellar medium (ISM) created by the moving neutron star and its wind), confined by ram pressure. X-ray emission is due to synchrotron emission from the wind particles accelerated at the termination shock, which is typically seen (if angular resolution permits) as the brightest portion of the extended structure (see, e.g., Kargaltsev & Pavlov 2008a), as predicted by MHD simulations (Bucciantini 2002; Van der Swaluw 2003; Bucciantini et al. 2005). In order to produce a population of high-energy particles, a highly energetic pulsar is required ($\dot{E}$ greater than ∼10³⁴; see, e.g., Kargaltsev & Pavlov 2008b). For classical synchrotron nebulae, we expect a relatively bright diffuse emission surrounding the pulsar, where the emission from the wind termination shock is brightest, as observed in all the other known cases (see, e.g., Gaensler et al. 2004; McGowan et al. 2006; Kargaltsev & Pavlov 2008b). The brightness is only slightly dependent on the ISM density ($\propto {n}_{\mathrm{ISM}}^{1/2}$); thus, we expect a relatively uniform brightness profile. One of the newly discovered pulsars, PSR J0357+3205 (Abdo et al. 2009), revealed a bright tail that cannot be interpreted in the framework of synchrotron PWNe, owing to the low energetics of the leading pulsar and to a peculiar shape (De Luca et al. 2011). Alternatively, a new model for pulsar nebulae was proposed, involving bremsstrahlung thermal emission from the ISM, shocked by the supersonically moving pulsar, accounting for all the characteristics of the nebula (Marelli et al. 2013). The cooldown time for this emission is orders of magnitude longer than for the synchrotron model, thus implying a lack of spectral variation along the tail. Pulsar bremsstrahlung nebula (PBN) emission is much more dependent on the ISM density ($\propto {n}_{\mathrm{ISM}}^{2}$), thus pointing to a clumpy brightness profile, with the brightness following and enhancing the medium distribution map. For this type of nebula we expect that most of the energy in the preshock flow is carried by the ions, while the electron temperature is responsible for the X-ray emission: we do not expect X-ray emission near the pulsar, owing to the Coulomb heating time of electrons. Both PWNe and PBNe obviously require the proper motion of the pulsar to be aligned with the main-axis nebular emission, in the case of a cometary shape. Recently, a peculiar nebular feature was found around the pulsar IGR J11014–6103 (Pavan et al. 2011). A deep Chandra observation revealed a helical structure protruding from the point source, possibly associated with a jet-like emission. A smaller trail was detected aligned with the direction of the putative, very large (>10³ km s⁻¹) proper motion, and it was associated with a classical pulsar wind nebular emission. This is reminiscent of the jet of the Guitar Nebula (PSR B2224+65, Hui & Becker 2007). This feature is usually explained in terms of a collimated ballistic jet, even if more complex explanations have been developed (Bandiera 2008). The helical shape and the collimated structure are the distinctive traits of this type of pulsar jet nebula. Another fundamental characteristic is the misalignment with the direction of proper motion.

Among the extremes of the population of Fermi pulsars, PSR J2055+2539 (J2055 hereafter; Saz Parkinson et al. 2010) is one of the best candidates. It is one of the 100 brightest γ-ray sources, with a γ-ray flux of (5.5 ± 0.1) × 10⁻⁷ erg cm⁻² s⁻¹ (Acero et al. 2015). It is located far from the Galactic plane, at a latitude ∼−13°, making the search for counterparts at other wavelengths easier. A blind search using Fermi data led to the unambiguous detection of the timing signature of a pulsar with P ∼ 320 ms and $\dot{P}\;\sim$ 4.1 $\times {10}^{-14}$ s s⁻¹ (Abdo et al. 2013). J2055 is among the less energetic and oldest pulsars in the Fermi zoo, with a spin-down energy of $\dot{E}$ = 5.0 × 10³³ erg s⁻¹ and a characteristic age ${\tau }_{c}$ = 1.2 Myr. The γ-ray "pseudo-distance" relation (Saz Parkinson et al. 2010), based on an empirical relation among γ-ray pulsar luminosities and their spin-down energies, yields an estimated distance of 600 pc. The second Fermi pulsar catalog (Abdo et al. 2013) reports J2055 to be one of the very few γ-ray pulsars with a strong detection of off-pulse emission, characterized by a spectral cutoff: this emission probably originates from its magnetosphere, opening a new conundrum for magnetospheric emission models. In fact, for young pulsars this is a serious challenge to outer-magnetosphere models radiating only above the null-charge surface; such weakly pulsed emission should be rare, being expected only for nearly aligned pulsars.

Owing to its peculiar characteristics, after two short XMM-Newton exploratory observations in 2013 we obtained a long, uninterrupted XMM-Newton observation to study the X-ray counterpart of J2055 and to search for X-ray pulsations. Section 2 reports the characteristics of the observations and the analysis methods. In Section 3 we report the analysis of the serendipitous sources in the field. Section 4 describes the spectral and timing characteristics of the pulsar in the X-ray band. This campaign revealed the presence of spectacular nebular features, described in Section 5. Pulsar and nebular modeling are treated in Sections 6 and 7. In the appendices, we describe the details of the field analysis (Appendix A) and of the new script we developed for our XMM-Newton analysis (Appendix B).

Our deep XMM-Newton observation of J2055, lasting 136.2 ks, was performed on 2013 May 1 (ObsID 0724090101). The PN camera (Struder et al. 2001) of the EPIC instrument was operating in Large Window mode, with a time resolution of 47.7 ms on a 27' × 13' field of view (FOV). The high time resolution, combined with the large FOV, allows for both the timing analysis of the J2055 pulsar and the spatial analysis of the nebular structures. The Metal Oxide Semiconductor (MOS) detectors (Turner et al. 2001) were set in full-frame mode (2.6 s time resolution on a 15' radius FOV). The thin optical filter was used for both PN and MOSs. We also analyzed XMM-Newton observations 0605470401 and 0605470901, taken on 2009 October 26 and on 2010 April 21, and lasting 24.5 and 17.9 ks, respectively. The PN camera was operating in full-frame mode, with a time resolution (73.4 ms). A timing analysis for this pulsar would be possible even using full-frame PN data, but it would result in a lower resolution of the light curve. We therefore used data from all the observations for the following analysis, with the exception of the timing analysis, for which we considered only PN data of the longest observation.

We used the XMM-Newton Science Analysis Software (SAS) v13.5. For the standard analysis, a tool included in the SAS epproc usually randomizes the arrival time of each event within the time resolution. Here we chose not to randomize the arrival times in order to achieve a better timing resolution for our pulsar. A discussion on the implication of this choice is reported in Marelli et al. (2014b). We performed a standard analysis of high particle background with the SAS tool bkgoptrate (also used for the compilation of the 3XMM source catalog⁶ ). This tool searches for the point at which the maximum signal-to-noise ratio (S/N) is achieved for the given background time series: time bins above this threshold are excluded. We then cross-checked the results using an independent method, following De Luca et al. (2005). Both the analyses revealed a high contamination from flares during almost the entire observation 0605470901, which was therefore excluded from subsequent analysis. Moderate contamination was detected in the other observations, lowering the good time interval to 132 ks. Following the prescription from the 3XMM source catalog, for the PN camera we selected 0–4 pattern events in the 0.5–10 keV energy range and 0 pattern events in the 0.3–0.5 keV energy range. We applied standard flags to remove bright columns. We selected 0–12 pattern events for the MOS detectors in the 0.3–10 keV energy range, applying a similar reduction for bright columns. Finally, we excluded out-of-FOV events and residual contaminants. Figure 1 shows images of the counts satisfying our event selection criteria from all the data sets. For the timing analysis, we used the SAS tool barycen to barycenter the PN events at the best-fit XMM-Newton pulsar position (see Section 3). For each spectrum, we generated ad hoc response matrices and effective area files using the SAS tools rmfgen and arfgen.

Zoom In Zoom Out Reset image size

We checked for other X-ray observations of the source and found a Suzaku (Mitsuda et al. 2007) observation (ObsID 405015010) taken on 2010 October 29 lasting 31.1 ks. The faintness of the X-ray counterpart and the crowding of the field (see Appendix A), coupled with the low spatial resolution of Suzaku (the half-power diameter of the XIS instrument, on board Suzaku, is ∼2'), make this observation useless for our purposes. It was thus excluded from our analysis.

Figure 2 shows the 0.3–10 keV XMM-Newton FOV. Source detection using maximum likelihood fitting was done both simultaneously and on each of the EPIC-PN, MOS1, and MOS2 data sets, using the SAS tool edetect_chain. The resulting detections were then confirmed using the CIAO tool wavdetect, following prescriptions from the CIAO online cookbook.⁷

Zoom In Zoom Out Reset image size

The best position of the pulsar X-ray counterpart is 20^h55^m4889+25°39'578 (03 + 15 1σ statistical plus systematic errors). From the most recent compilation of the LAT γ-ray pulsar timing models,⁸ the best γ-ray position of the pulsar is 20^h55^m4894+25°39'591 (05 1σ statistical plus systematic error), fully compatible with the X-ray position. According to the log N–log S distribution of the Chandra galactic sources (Novara et al. 2009) at mid-Galactic latitudes, we can estimate the probability of a chance detection of an X-ray source within the LAT γ-ray timing error box and with X-ray flux comparable to or greater than the X-ray counterpart ($\sim {10}^{-14}$ erg cm⁻² s⁻¹) to be less than 10⁻⁵. Thus, based on positional coincidence, we consider our identification secure.

In Figure 2 we marked the brightest sources in the PN FOV. A bright nebular source is apparent in the field. This object, possibly associated with a galaxy cluster, will be described in F. Gastaldello et al. (2016, in preparation) and not further described in this paper. We chose all the sources with an absorbed flux greater than 5 $\times {10}^{-15}$ erg cm⁻² s⁻¹, with a 5σ detection and with at least 150 XMM-Newton net, total counts. We also considered all the sources detected in the nebular regions. The study of the line-of-sight absorption of such sources could allow us to constrain the pulsar distance. Indeed, after selecting candidate active galactic nuclei (AGNs) in the FOV, it is possible to measure, from their spectra, the total Galactic column density in the direction of J2055. Next, using the pulsar column density, we can get an estimate of its distance with respect to the edge of the Galaxy. Such an estimate could be refined if bright X-ray-emitting and optical-emitting stars (with known distances) were also present (see, e.g., Marelli et al. 2013, 2014a, 2014b). The spatial resolution of this measurement is much higher than the one from the 21 cm ${{\rm{H}}}_{{\rm{I}}}$ sky survey of Kalberla et al. (2005). The resolution of the survey is 06, larger than the entire XMM-Newton FOV: the presence of local molecular clouds could affect the ${{\rm{H}}}_{{\rm{I}}}$ result even by orders of magnitude (see, e.g., the case of the PSR J1813–1246 field; Marelli et al. 2014b).

Based on the study of the spectra of serendipitous sources and of associated optical counterparts, we classified serendipitous sources as AGNs or candidate stars (see Appendix A). Our exercise allowed us to identify 24 very likely AGNs and 13 stars from a total of 58 X-ray sources. In particular, one of the brightest sources (source 1 in Table 1) was originally detected as a single source: only through the application of a new script were we able to determine that it consisted of the superimposition of multiple sources (see Appendix B for details).

Table 1.  Analysis of the Brightest Serendipitous Sources in the XMM-Newton FOV

Source	J2000 Coordinates	Spectrum	nH	$\mathrm{log}\left(\frac{{f}_{{\rm{X}}}}{{f}_{\mathrm{opt}}}\right)$ a	ID
			$({10}^{21}\;{\mathrm{cm}}^{-2})$
1	313.96143 25.679886	apec+pow/2apec	${4.77}_{-1.21}^{+0.78}$/4.00${}_{-0.27}^{+0.53}$	>0.566 + 0.686/>1.12(V*)	star + agn
2	313.94955 25.687331	2apec	${0.981}_{-0.103}^{+0.118}$	−1.58(V)	star
3	313.94435 25.694125	pow	${65.5}_{-11.0}^{+12.8}$	?	agnb
4	313.93596 25.706747	pow	${1.23}_{-0.48}^{+0.60}$	>−0.291(V*)	agn
6	313.94593 25.651206	pow/apec	${2.79}_{-0.36}^{+0.39}$/2.06 ± 0.25	−0.927/−0.503(V*)	star
7	313.95363 25.631689	pow/apec	${2.05}_{-0.48}^{+0.56}$/1.51${}_{-0.28}^{+0.34}$	>−0.255/>−0.884(V*)	?
8	313.95196 25.607475	pow	${2.28}_{-0.71}^{+0.86}$	>−0.546(V*)	agn
9	313.9496 25.590583	apec	${4.60}_{-0.69}^{+2.03}$	−3.86(V)	star
10	313.94093 25.547181	pow/apec	${2.06}_{-0.40}^{+0.44}$/1.52${}_{-0.25}^{+0.26}$	>−0.051/>0.332(V*)	?
13	313.96099 25.6253	pow	${0.135}_{-0.135}^{+0.298}$	−1.86(V)	?
14	313.97678 25.609422	pow/2apec	${1.80}_{-0.10}^{+0.11}$/3.65${}_{-0.81}^{+0.78}$	>0.607/>0.899(V*)	?
15	313.99079 25.606825	pow/2apec	${1.68}_{-0.31}^{+0.35}$/2.63${}_{-1.03}^{+3.29}$	−0.387/−0.612(B)	star
16	313.97448 25.586575	apec	${1.18}_{-0.53}^{+0.81}$	−3.29(V)	star
17	314.00763 25.609944	pow	${1.45}_{-0.27}^{+0.31}$	>0.127(V*)	agn
18	313.96974 25.538703	pow	${1.57}_{-0.59}^{+0.81}$	>−0.367(V*)	agn
19	313.98956 25.539306	pow	${1.33}_{-0.34}^{+0.42}$	>−0.129(V*)	agn
20	313.92683 25.654367	pow	${3.44}_{-0.50}^{+0.56}$	>−0.273(V*)	agn
21	313.92368 25.690575	pow	${2.63}_{-0.36}^{+0.40}$	>−0.108(V*)	agn
22	313.98105 25.698697	pow	${0.366}_{-0.143}^{+0.207}$	>0.297(V*)	agn
23	313.98844 25.702833	pow/apec	${5.26}_{-1.24}^{+1.56}$/4.04${}_{-0.97}^{+1.50}$	>−0.506/>−0.0326(V*)	?
24	313.95607 25.707142	pow/apec	${3.67}_{-1.43}^{+2.39}$/2.41${}_{-0.96}^{+1.68}$	>−0.695/>−0.166(V*)	?
25	313.96015 25.712533	pow	${15.6}_{-5.42}^{+6.75}$	>−2.47(V*)	agn
26	313.90219 25.666081	pow	${0.164}_{-0.08}^{+0.09}$	>−0.0620(V*)	agn
27	313.86733 25.666336	pow	${1.56}_{-1.51}^{+4.25}$	>−0.310(V*)	agn
28	313.84279 25.631844	pow	${0.816}_{-0.339}^{+0.527}$	>−0.321(V*)	agn
29	313.80258 25.653275	pow	${1.49}_{-0.33}^{+0.40}$	−1.35(V*)	?
30	313.79122 25.650547	pow	${2.08}_{-0.31}^{+0.34}$	>0.160(V*)	agn
31	313.74973 25.656481	pow/apec	${8.12}_{-2.25}^{+2.93}$/7.48${}_{-2.35}^{+2.57}$	>−0.929/>−0.597(V*)	?
32	313.88097 25.580411	pow/apec	${3.32}_{-1.41}^{+1.89}$/2.05${}_{-0.78}^{+1.14}$	>−0.743/>−0.368(V*)	?
33	313.84338 25.568675	?	?	?	?
34	313.77725 25.530489	pow	5.50	−1.56(V)	?
35	313.89859 25.535194	pow	${2.33}_{-0.43}^{+0.51}$	>−0.0162(V*)	agn
36	313.91509 25.61615	pow/2apec	${1.71}_{-0.36}^{+0.42}$/1.32${}_{-0.34}^{+0.77}$	>−0.192/>0.187(V*)	?
37	314.00874 25.573114	pow/apec	${4.17}_{-2.00}^{+3.23}$/4.13${}_{-1.80}^{+3.78}$	−1.66/−1.44(V*)	star
39	314.05828 25.556017	pow/apec	${5.53}_{-1.38}^{+2.00}$/4.80${}_{-1.41}^{+2.24}$	>−0.540/>−0.147(V*)	?
40	314.03453 25.595597	pow	${1.32}_{-0.19}^{+0.22}$	>0.176(V*)	agn
41	314.02360 25.612442	pow	<0.629	>−0.277(V*)	agn
42	314.04116 25.624814	pow	${2.15}_{-0.90}^{+1.22}$	>−0.653(V*)	agn
43	314.01680 25.684672	pow	${8.10}_{-2.51}^{+3.33}$	>−1.35(V*)	agn
44	314.03226 25.693653	apec	${1.42}_{-0.53}^{+0.72}$	−2.50(V*)	star
45	314.04296 25.708472	pow/apec	${7.75}_{-2.02}^{+2.61}$/5.46${}_{-1.24}^{+1.88}$	>−0.938/>−0.465(V*)	?
46	314.06314 25.700192	pow/apec	${1.32}_{-0.12}^{+0.13}$/0.977${}_{-0.084}^{+0.079}$	−0.288/0.0897(V*)	star
47	314.06771 25.673589	apec	${4.64}_{-1.48}^{+1.50}$	−3.16(V)	star
48	314.09986 25.658803	pow/apec	${1.74}_{-0.78}^{+1.03}$/1.12${}_{-0.63}^{+0.40}$	>−0.435/>−0.146(V*)	?
49	314.10914 25.736881	pow	${1.63}_{-0.38}^{+0.45}$	>0.0333(V*)	agn
50	314.03943 25.747156	pow	${1.87}_{-0.10}^{+0.11}$	0.423(B)	agn
51	314.01804 25.758858	pow	${3.77}_{-0.79}^{+0.97}$	−0.0679(V*)	agn
52	313.95291 25.776819	pow	${2.35}_{-0.14}^{+0.15}$	0.403(B)	agn
53	313.90073 25.751294	apec	${1.09}_{-0.36}^{+0.47}$	−2.12(V)	star
54	313.93571 25.858336	pow	${2.23}_{-0.41}^{+0.46}$	0.0676(V*)	agn
55	313.77492 25.716167	apec	${1.49}_{-0.39}^{+0.49}$	−3.74(V)	star
56	313.95314 25.445281	pow/apec	<1.00/<0.420	>0.118/>0.448(V*)	?
57	313.96553 25.421706	pow	${1.34}_{-0.94}^{+1.79}$	>0.000970(V*)	agn
58	313.97756 25.717731	pow	${0.249}_{-0.228}^{+0.596}$	−1.12(V*)	?
59	313.99449 25.719272	pow/apec	${6.10}_{-1.24}^{+1.50}$/4.73${}_{-1.12}^{+2.59}$	>−0.722/>−0.345(V*)	agnb
60	313.90637 25.638078	pow/apec	${2.17}_{-0.40}^{+0.46}$/1.38${}_{-0.29}^{+0.33}$	−1.43/−1.93(V)	star
61	313.90309 25.630814	pow	${1.93}_{-0.30}^{+0.38}$	−0.733(B)	?

Notes. The table shows the best source position, the statistically acceptable spectral model fit(s), the fitted column density(s), the X-ray-to-optical flux ratio, and the resulting identification. Typical XMM-Newton positional errors are 3''. We reported the 1σ error on column density.

^aOptical-to-X-ray flux ratio for the different fitting models. Where present, we used the V-band optical flux. If not detected, we used the B-band optical flux. In case of an R-band detection, we obtained the V-band flux (see Appendix A). ^bIdentification of these sources is based on optical and hardness ratio tests.

Download table as:  ASCIITypeset images: 1 2

The simultaneous fit of AGN spectra, with the column density linked, results in a relatively low value of column density: (2.07 ± 0.08) $\times \;{10}^{21}$ cm⁻² (1σ error), as expected in a midlatitude region. This result is higher than the value of 1.2 × 10²¹ cm⁻² obtained from the 21 cm ${{\rm{H}}}_{{\rm{I}}}$ sky survey of Kalberla et al. (2005).

In order to optimize the pulsar extraction radius, we evaluated the brightness profile of the putative counterpart to maximize both the S/N and the number of counts. We found the best extraction radius to be $15^{\prime\prime}$, containing a total of 1020 and 636 counts in the 0.3–10 keV energy range in the PN and the two MOSs detectors, respectively, with a background contribution of ∼25%. We generated ad hoc response matrix and effective area files using the SAS tools rmfgen and arfgen. Owing to the relatively low statistics and to optimize the fit, for the spectral analysis we added all the data from the MOS cameras using the HEAsoft tools mathpha, addarf, and addrmf to sum spectra, effective area files, and response matrix, respectively. PN spectra were analyzed separately owing to the different operation modes (large window and full frame) requiring different response matrices. The background region was extracted from a nearby source-free region common for every camera and on the same CCD of the source. For the pulsar emission, we tried a simple power-law (pow), a simple blackbody (bb), a simple magnetized neutron star atmosphere (nsa), the combination of a power-law and a blackbody (bb), and the combination of a power-law and a magnetized neutron star atmosphere model (nsa in XSPEC—assuming a neutron star with a radius of 13 km, mass of 1.4 M_⊙, and surface magnetic field of 10¹³ G). The simple power-law description gives an acceptable fit (reduced chi square ${\chi }_{\mathrm{red}}^{2}=1.05,$ 61 degrees of freedom (dof)) with a column density of N_H = (2.18 ± 0.26) × ${10}^{21}$ cm⁻² and a photon index of Γ = 2.36 ± 0.14 (1σ confidence errors). The best-fit unabsorbed 0.3–10 keV flux is (3.43 ± 0.27) × 10⁻¹⁴ erg cm⁻² s⁻¹. This flux is in agreement with the one reported in Abdo et al. (2013) of ${4.3}_{-2.8}^{+1.2}\;\times$ 10⁻¹⁴ erg cm⁻² s⁻¹. Figure 3 shows the spectra with the best-fit power-law model. Neither the single-component blackbody nor the nsa model gives an acceptable fit (${\chi }_{\mathrm{red}}^{2}$ = 2.05 ${\chi }_{\mathrm{red}}^{2}$ = 1.81, respectively, 61 dof), and both were therefore rejected. Although the combined power-law plus thermal models give acceptable fits (both the bb and nsa give ${\chi }_{\mathrm{red}}^{2}$ = 0.99, 59 dof), the thermal normalization is compatible with zero at 2σ. Moreover, an F-test (Bevington 1969) performed comparing the simple power law with the composite spectra does not point to a significant improvement by adding a thermal component. We thus consider the single-component nonthermal model to be the best for our source.

Zoom In Zoom Out Reset image size

To search for X-ray pulsations, we used the most recent Fermi ephemeris.⁹ As reported in Section 2, we used only barycentered PN events from the 2013 XMM-Newton observation. We decided to apply a photon-weighting technique similar to the one used for Fermi-LAT (using gtsrcprob; Kerr 2011), assigning to each photon a probability of coming from the pulsar. This allows for a better background rejection and improves the sensitivity to pulsations, as discussed in Marelli et al. (2014b). To this end, we further developed the Python tool reported in Marelli et al. (2014b) (see Appendix B). Considering only photons with a probability greater than 0.01 of coming from the pulsar, we used a weighted Markov Chain Monte Carlo algorithm (MCMC; see, e.g., Wang et al. 2013) to test and refine the γ-ray ephemeris. Owing to the shortness of the X-ray observation, we fixed the derivatives of frequency to ones of the γ-ray ephemeris. We obtain a best-fitted frequency F = 3.129286 ± 0.000003 s⁻¹, in agreement with the γ-ray one. With six harmonics, the resulting H-value of 59.09 gives a tail probability of 1.25 × 10⁻¹¹. Even taking into account the ∼50 trials, we found the X-ray pulsation with a confidence greater than 6σ (see Figure 4 for the phased light curve), thus further confirming the identification of our X-ray candidate with the γ-ray pulsar. Unfortunately, the low statistics hamper a detailed energy-dependent timing analysis as was done in the case of the CTA-1 pulsar (Caraveo et al. 2010). This could have allowed us to unambiguously accept or exclude a possible low-energy (e.g., thermal) component of the spectral model.

The most significant feature in the XMM-Newton FOV (Figure 2) is a thin, straight, extended emission passing through half of the detectors north to south and protruding from the position of J2055. This feature is detected in each data set, so that we exclude an instrumental origin. Much less obvious is another elongated extended emission protruding from J2055 that is visible southeast to northwest. Unfortunately, this fainter secondary feature is heavily contaminated by bright point sources, hampering a classical analysis (see, e.g., Marelli et al. 2013).

In order to better disentangle the nebular emission from the background and the numerous superimposed point-like sources, we developed a new Python-based tool to simulate and subtract their contribution from XMM-Newton data and images, based on their positions and spectra (see Appendix B). After running our tool, for each observation we created a new data set that excludes photons based on the probability, P, that they come from one of the considered sources (pulsar, serendipitous sources, and background): using a Poisson distribution, each photon has a probability P of being excluded. The remaining photons then come from either statistical fluctuations or other sources not considered, such as our extended features. We then based our analysis of extended sources on this reduced data set, subsequently confirming our results on the original data sets.

Zoom In Zoom Out Reset image size

5.1. Spatial Analysis

For each camera, we produced brightness profiles of each nebula along the main axis (AMA) of the feature and orthogonal to the main axis (OMA). At this point, we chose the main axis only through a visual evaluation, ending at the pulsar position. In order to maximize the statistical significance of each bin, we carefully evaluated the bin size and the width of extraction boxes, choosing a scale of 100'' × 50'' AMA and 1000'' × 10'' OMA for the main feature (shown in orange and magenta, respectively, in Figure 5) and 100'' × 50'' AMA and 400'' × 20'' OMA for the secondary feature (shown in yellow and white, respectively, in Figure 6). We extracted the number of counts in each box, summing the results of each data set. In order to take into account the spatial evolution of the exposures, as well as the lack of data from an entire camera at some points, we also summed the exposure maps, produced through the SAS tool eexpmap, and we built a graph of the exposure evolution AMA and OMA, using the same boxes as before. This graph was then normalized at the mean value of the exposure. Finally, we divided the counts profile by the exposure profile to obtain a normalized brightness profile, in normalized counts arcsecond⁻². Then, for each nebula we proceded iteratively by slightly changing the main axis in order to find the maximum of the OMA profile compatible with zero. Therefore, we will consider this as the main axis of each nebula and its error as the propagated error from the maximum of the OMA profile. The brightness profiles are shown for the main feature in Figures 7 and 8 and for the secondary feature in Figures 9 and 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

First, we tested the presence of the fainter nebula. To do this, we chose visually an ad hoc area containing most of the nebular photons. We considered different nearby regions with the same shape and on the same CCDs as the secondary nebula extraction area. After considering the different exposure, as well as instrumental effects, we obtain an expected number of counts from the background in the secondary nebula extraction region of ∼12,600. We estimated, conservatively, that ∼10% of the total counts are due to residual photons from point-like sources superimposed on the nebula, thus obtaining a total contribution of ∼800 counts. We detected a total of 15,130 counts from the secondary nebular region (after taking into account all the instrumental effects). The tail probability of obtaining 15,130 counts from a Poisson distribution peaked at 13,400 confirms the presence of this feature at more than 7σ. A similar calculation performed on the raw data, with 95% of the counts from point sources removed using circular regions, gives as well a nebular detection at more than 5σ. We considered the possibility that such an excess comes from point-like sources not detected by the source detection algorithms described in Section 3. These cannot realistically account for such a high number of additional counts: by taking into account the upper limit on the detectable flux, we expect no more than ∼100 total counts from each of them, thus requiring a large number of unresolved sources.

We obtain the main axis of the main nebula to be at −01 ± 02 with respect to the north–south axis, counterclockwise. In the same reference, the secondary nebula axis is at 1627 ± 07 (1σ errors).

As shown in Figure 7, AMA the main nebula is characterized by a slow increase in flux until ∼9', followed by a rapid drop until it becomes undetectable at ∼12'. Orthogonal to the main axis (Figure 8), the nebular profile is clearly peaked, following a Lorentzian-like distribution that fades symmetrically to the axis, ∼20'' from it. We note that the OMA distribution is, at the limit, compatible with the theoretical XMM-Newton PSF distribution (${\chi }_{\mathrm{red}}^{2}$ = 2.8, 11 dof): this points to a thin (no more than ∼20'' thick), straight nebula or a thinner, curved nebula. To detect a possible main feature curvature with respect to the main axis, we extracted counts from boxes of 170'' × 125'' AMA. Each box was than divided into 15'' × 125'' sub-boxes, following the same procedure as before for the exposure, and then we fitted each box with a Gaussian. All the fits are statistically acceptable, and the values (and errors) of the maximum position are reported in Figure 7. The resulting curve is not compatible with a constant model (${\chi }_{\mathrm{red}}^{2}$ = 4.2, 6 dof): the main nebula shows a curvature in the clockwise direction at ∼3' and then another in the counterclockwise direction at ∼9'.

Along the main axis, the secondary nebula is characterized by an almost flat flux until ∼150'', then slowly increasing until ∼250'', where it suddenly fades. Orthogonal to the main axis, it is peaked following a Gaussian-like distribution with a sigma of 27'' ± 5'' (1σ error). It does not follow the theoretical XMM-Newton PSF distribution (${\chi }_{\mathrm{red}}^{2}$ = 17.3, 6 dof): it is thicker (at least 50'' thick) than the main nebula. The low statistics prevent us from carrying out a more detailed spatial analysis.

Brightness profiles of the tails in the soft (0.3–1.5 keV) and hard (1.5–6 keV) energy bands show no statistically significant differences from the total profiles.

5.2. Spectral Analysis

For the spectral analysis, we evaluated the best nebular extraction regions from the results of the spatial analysis. For each nebula, we chose extraction boxes that ensure that the nebula is seen with more than a 3σ significance. This results in a box of 700'' × 50'' for the main feature and a box of 250'' × 80'' for the secondary feature, oriented as their main axes. These boxes were also divided into two equal parts AMA in order to study the nebular spectral evolution. For spectral analysis, we only took into account cameras that cover more than 50% of the considered region. The normalizations of nebular spectra were left free to vary for each camera, owing to the different FOV and coverage of the instruments. We extracted background counts from regions with the same shape and angle, passing through the same CCDs as the source regions. We chose these ad hoc regions to be as near as possible to the source regions, but excluding the nebular features and the extended feature west of the main nebula. In order to be conservative, we removed from the region 5''-radius circles centered at the position of point-like sources, thus excluding possible residuals. We followed the same procedure as in Section 4 to produce spectra, but we chose not to add them owing to the higher statistics. Spectra were grouped in order to have at least 150 counts per bin: owing to the high background, this allows us to have at least 25 source counts per bin in each camera.

A power-law model fits well both nebulae. The best-fit absorption columns of the nebulae and pulsar are compatible. Then, under the hypothesis of association between the nebulae and the pulsar (see Section 6 for more details), we fitted all the spectra together, linking the column density. In this way, we find the column density of the system to be (2.22 ± 0.17) × 10²¹ cm⁻² (90% confidence error). We note that a thermal bremsstrahlung model fits well both the nebulae, with column densities again in agreement with the one of the pulsar. We found no sign of spectral evolution with distance from the pulsar for both the nebulae down to a variation of 0.1 and 0.3 (90% confidence) in the photon index of the main and secondary nebula, respectively. See Table 2 for details on the fits.

Table 2.  Results from Spectral Analysis of Pulsar and Nebulae

Source	NH	Γ/kT	Fluxa	${\chi }_{\mathrm{red}}^{2}$/dof
	(1021 cm−2)	(−/keV)	(10−14 erg cm−2 s−1)
Total	2.22 ± 0.17	...	...	...
psr	2.18 ± 0.26	2.36 ± 0.14	3.43 ± 0.27	1.05/61
pwn1	2.25 ± 0.23	1.82 ± 0.08	28.7 ± 1.9	1.19/249
pwn1-br	1.68 ± 0.16	${5.41}_{-0.73}^{+0.92}$	23.5 ± 1.8	1.15/249
pwn2	${2.11}_{-0.59}^{+0.73}$	1.62 ± 0.21	5.45 ± 1.22	1.14/139
pwn2-br	${1.60}_{-0.40}^{+0.49}$	10.3 ± 4.4	4.85 ± 0.84	1.15/139
pwn1a	...	1.75 ± 0.12	...b	1.11/786
pwn1b	...	1.82 ± 0.11	...	...
pwn2a	...	1.50 ± 0.33	...	...
pwn2b	...	1.77 ± 0.20	...	...

Notes. All the errors are at 1σ confidence level. Here we report the column density obtained by fitting the pulsar and the two nebular spectra simultaneously (total). We report the pulsar spectrum. We report the spectra of each of the nebulae (the main, 1, and the secondary, 2). Lastly, we report the contemporaneous spectra of the two parts (nearer to pulsar, a, and farthest from pulsar, b) of the two nebulae, with the column density linked.

^aUnabsorbed 0.3–10 keV flux. Here we report the flux of the instrument with the more complete view of the nebula (respectively MOS 2 of the second observation and PN of the second observation for the main and secondary nebula). ^bThe flux of the two parts of the nebulae is useless in order to describe them. For a complete spatial analysis see Section 5.1.

Download table as:  ASCIITypeset image

6.1. The Pulsar

6.2. The Nebulae

We analyzed the XMM-Newton observations of the radio-quiet γ-ray PSR J2055+2539. We confirmed the X-ray counterpart association of the pulsar. The X-ray emission is nonthermal, with an X-ray efficiency of ${\nu }_{{\rm{x}}}$ = L_x/$\dot{E}\;\sim$ 3 × 10⁻⁴ (for a distance of 600 pc), typical of X-ray pulsars. The X-to-γ-ray luminosity ratio is L_x/${L}_{\gamma }\;\sim$ 1600, typical of radio-quiet pulsars. The lack of cooling thermal emission is due to the old age of the pulsar, while the lack of hot spot thermal emission can be due to its low energetics. We detected nonthermal X-ray pulsations, following a sinusoidal profile with a pulsed fraction of ∼25%. Taking into account considerations on the γ-ray efficiency of the pulsar and on its X-ray spectrum, we can infer a pulsar distance ranging from 450 to 750 pc.

We developed a new script to simulate the spatial distribution of XMM-Newton counts from input point-like and extended sources. Using this, we found and analyzed two different nebular features associated with J2055. Both nebulae have an elongated profile, protruding from the pulsar and almost symmetrical to a main axis. The main axes of the two nebulae are separated by an angle of ∼163°. The main, brighter feature is 12' long, <20'' thick (corresponding to 2.1 × 0.05 pc at a 600 pc distance), characterized by an asymmetry with respect to the main axis evolving with the distance from the pulsar, possibly forming a helical pattern. The secondary feature is 250'' × 30'' (which corresponds to 0.7 × 0.09 pc at a 600 pc distance). Both nebulae present an almost flat brightness profile, with a sudden decrease at the end. The nebulae can be fitted by either a power-law model or a thermal bremsstrahlung model. Their efficiencies (L_x/$\dot{E}$) are 2.4 × 10⁻³ and 3.0 × 10⁻⁴ for the main and the secondary nebula, respectively. A plausible interpretation of the brighter nebula is in terms of a collimated ballistic jet, mainly owing to its thickness and shape. The secondary nebula is most probably a classical synchrotron-emitting tail. Future Chandra observations will allow us to better study the shape of the nebulae and to find the proper motion of the pulsar, thus allowing us to test our interpretations. Finally, we found and studied a bright, serendipitous galaxy cluster with a central AGN.

This work was supported by the ASI-INAF contract I/037/12/0, art.22 L.240/2010 for the project "Calibrazione ed Analisi del satellite NuSTAR."

Facility: XMM - Newton X-Ray Multimirror Mission satellite.

Based on the study of spectra and possible optical counterparts, we can classify serendipitous sources as AGNs or candidate stars, allowing us to constrain the pulsar distance. Indeed, after selecting candidate AGNs in the FOV, it is possible to measure from their spectra the total Galactic column density in the direction of J2055.

For this analysis, we decided to use only brighter sources, those from which it is possible to extract and fit a spectrum: we therefore require that they have an absorbed flux greater than 5 $\times \;{10}^{-15}$ erg cm⁻² s⁻¹, with a 5σ detection and with at least 150 XMM-Newton net, total counts. We also considered all the sources detected in the nebular regions. As described in Section 3, we found 57 serendipitous bright X-ray sources in the FOV. In order to classify them, we performed a spectral analysis, as done in Marelli et al. (2014b). All the sources, with the exception of source 1, are compatible with a point-like distribution. As described in Appendix A, a deeper analysis of source 1 revealed it as the superimposition of two or three different sources. Extended sources are discussed elsewhere. As a first method of classification, we used their spectra. Based on Novara et al. (2009), at mid-Galactic latitude we expect that ∼75% of serendipitous sources with the considered flux are AGNs and the remaining are stars. X-ray emission from AGNs is best described by a power law, while X-ray emission from stellar coronae by single or double apec (emission spectrum from collisionally ionized diffuse gas). For sources located inside extended emissions, we extrapolated the spectrum and flux of the extended source and froze it into the spectral fit. Out of these 58 sources, 30 can be realistically fitted only be a power law and 5 only by apec(s); for the remaining 23 objects, this type of analysis is not conclusive (see Table 1).

Since the count statistics of some of the selected X-ray sources are too low to discriminate the spectral model, we performed a qualitative spectral analysis using the count rate (CR). The CR is measured in three energy ranges: soft (0.3–1 keV), medium (1–2 keV), and hard (2–10 keV). We then computed the two different hardness ratios:

$$\begin{eqnarray*}\mathrm{HR}12 & = & [\mathrm{CR}(1-2)-\mathrm{CR}(0.3-1)]/[\mathrm{CR}(1-2)\\ & & +\mathrm{CR}(0.3-1)]\\ \mathrm{HR}23 & = & [\mathrm{CR}(2-10)-\mathrm{CR}(1-2)]/[\mathrm{CR}(2-10)\\ & & +\mathrm{CR}(1-2)].\end{eqnarray*}$$

Large/small HR12 values point to a large/small absorption; a large/small HR23 indicates a hard/soft spectrum. The values of the hardness ratios are compared to a matrix of simulated values for apec (kT ranging from 0.5 to 5.5 keV) and power-law (Γ ranging from 1.5 to 2.5) models. Each spectral model is computed using the average ISM absorption given by Dickey & Lockman (1990). We extracted and simulated CRs from the PN camera of the second observation and the sum of MOS cameras, and then we combined them. These values are plotted in Figure 11 together with a sample of the measured HR ratios. From the position of the latter with respect to the former, an interpolation is possible over the two models and the values of the model parameters. Although this method proved best for the case of the highly absorbed serendipitous sources in the field of PSR J1813–1246 (Marelli et al. 2014b), for J2055 field this method did not improve the spectral results already obtained for most of the sources. Only in the case of source 2 do we have an indication of a star-like spectrum and for sources 3, 43, and 59 of AGN-like.

Zoom In Zoom Out Reset image size

Another common way to confirm X-ray classification of sources is based on multiwavelength analysis: the X-ray-to-optical flux ratio is a good indicator of the nature of X-ray emitters. According to La Palombara et al. (2006), AGNs have typical logarithms of X-ray-to-optical flux ratios higher than −0.15 (+0.2), while stars lower than +1.0 (+0.65) in the V (B) band. Inside each X-ray source error box, we looked for association with optical sources from the NOMAD catalog (Zacharias et al. 2005), considering the V-band magnitude as reference or B-band where necessary. Alternatively, we used the R-band magnitude: V-band magnitude was then evaluated from R using the following method. We considered the best-fitting X-ray models to estimate the column densities for each source. In case of multiple fitting models, we repeated the following for each model. We de-absorbed R magnitudes by using coefficients from Predehl & Schmitt (1995) and Cardelli et al. (1989), obtaining the R-band absolute magnitude. Then, assuming a thermal spectrum and/or a power-law spectrum, based on the X-ray model, we used the online NICMOS tool¹⁰ to evaluate the V-band absolute magnitude. Few of our AGN-like objects have optical counterparts inside their error box that could be due to spurious matches. In order to estimate the number of spurious matches, we used the relation from Severgnini et al. (2005) and Novara et al. (2009). This yielded a probability of chance coincidence of 21%, which means that, at our limiting magnitudes, contamination effects cannot be ignored. Each of the objects revealed as star-like from spectral or HR analysis has an optical counterpart that agrees with the expected X-ray-to-optical flux ratio. Table 1 reports the associated optical counterparts and expected upper limits for AGNs.

Finally, we combined results from spectral and multiwavelength analysis to classify our serendipitous sources, with the restrictions reported above. All the sources with spectral model or HR compatible with a star coronal emission and with a star-like X-ray-to-optical flux ratio are considered as stars: there are 13 of these objects. All the sources with spectral model or HR compatible with an AGN emission and with an AGN-like X-ray-to-optical flux ratio are considered as stars: there are 24 of these objects. A total of 21 sources cannot be uniquely classified.

We developed a script to simulate the spatial distribution of photons from input point-like and extended sources. After this simulation, by setting the value of the variable "other," it is possible either (i) to search for new sources using the reduced data set without photons from considered sources (other = yes) or (ii) to extract photons from a single source to perform high-level spectral and timing analysis (other = no). We note that this type of approach is particularly well suited to combine the high spatial resolution of the Chandra-ACIS observatory with the high spectral and timing resolutions of XMM-Newton. The first satellite can provide the number and exact position of point-like sources, as well as the flux and shape of nebular features: these can be used as inputs of the script to perform high-level XMM-Newton analysis. In this paper, first we used the script to clean out the nebular features from point-like sources reported in Table 1; see Figure 12. Then, we used again the script to extract photons from the pulsar to perform the timing analysis.

Zoom In Zoom Out Reset image size

We used Python v2.7.5 language to build the script. This Python script recalls standard Python modules (os, sys, math, shutil, multiprocessing, subprocess), as well as modules dedicated to scientific analysis (xspec, numpy, pyfits). This script has a number of variable input parameters:

• 1.  
  the type of analysis: to search for new sources or to extract photons from a source ("other" variable);
• 2.  
  the energy range and energy resolution to be used, in eV;
• 3.  
  the spatial accuracy to be used for point-like PSF reconstruction; if one of the sources is extended, it should be set also to the extended PSF accuracy;
• 4.  
  the names and positions of input and output files.

The script takes a list of inputs from an ASCII file:

• 1.  
  the position of each source, in both ecliptic and physical coordinates;
• 2.  
  the spectral model and parameters of each source, readable by xspec module;
• 3.  
  the source extraction circle radius from which each spectrum of a point-like source is extracted;
• 4.  
  the spectrum of background, readable by the xspec module;
• 5.  
  the extraction circle radius of background spectrum.

Finally, as input it needs some FITS files:

• 1.  
  event file, filtered by the user for instrumental effects and energy;
• 2.  
  image file constructed from the event file;
• 3.  
  exposure file constructed from the event and image files;
• 4.  
  calibration file (for XMM-Newton, the CCF file);
• 5.  
  the spectrum, background spectrum, response matrix (rmf), and effective area (arf) file from each source and from the background;
• 6.  
  if extended sources are loaded, for each of them an image of the theoretical spatial shape of the nebula (before the PSF application), possibly with different normalizations of each pixel.

As output, it produces the following:

• 1.  
  for each input source, an ASCII file containing all the commands ran by the script;
• 2.  
  for each input source and the background, a column is added to the fits file: this shows the probability (0 < P < 1) of the photon coming from that source;
• 3.  
  for each input source, the image of the simulation of that source (see Figure 12);
• 4.  
  an image containing all the simulated sources, an image with residuals, and an image containing the deviation of the simulated map from the counts map, in sigmas; owing to the high computing time needed, each source is treated separately on a different CPU and the results are merged only at the end, to compute probabilities and create images.

The first part of the script produces the mean CR over the entire observation of each source, in each energy band. Here the spectrum, background spectrum, rmf, and arf are loaded using the pyxspec library, and the input spectral model is assigned (but not fitted). Based on the spectral model interpolated with rmf and arf, the script computes the CR. We note that these values are corrected for PSF tails and exposure map.

The second part of the script constructs the PSF of each source, in each energy band. We used the last theoretical model for the XMM-Newton PSF, as reported in Read et al. (2011) and following calibration documents. The normalization of the PSF function is set in order to obtain the correct number of expected CRs, as calculated in the first part of the script. In case of an extended source, this is treated as a number of point-like sources with positions and normalizations following the input model image. We note that, in order to reduce the computing time, a different, lower-resolution, spatial grid can be set for extended sources. The use of this complicated, spatial, and energy-dependent PSF function rather than the simple King function previously adopted resulted in a significant improvement in the positional accuracy, as well as spurious source identification for automated XMM-Newton searches of serendipitous sources. This change in the PSF becomes more important for bright sources in the field, allowing for a better disentangling of photons from each source, as well as the detection of nearby faint sources. In the case of the J2055 field, the adoption of this newer PSF resulted in the clear identification of the two, or possibly three, different sources responsible for the emission from source 1 in Table 1.

The third part of the script computes the number of observed counts, from each source in each energy band, in each pixel of a spatial grid, with a step set by the user. It is based on the theoretical PSF function and values obtained in the second step. Finally, the obtained matrix is multiplied for the input exposure map to obtain the number of expected counts.

The fourth part of the script creates output images. The selected spatial grid is collapsed to agree with the input counts map resolution, and results from all the energy bands are summed. For each source, a simulated image is written, and they are all summed to write the total simulated sky. The script also creates a residual image, subtracting the simulated map (S) from the counts map (C), and a significance map defined for each i bin as SI_i = $\parallel$O_i − S${}_{i}\parallel$/$\sqrt{{O}_{i}}$, where O_i! = 0 and SI_i = 0 otherwise.

The fifth part of the script computes the probabilities. The script assigns to each photon in the FITS file a list of predicted counts values for each source, based on the pixel of the spatial grid where the photon falls. If "other" = yes, it also takes into account the observed number of counts in that pixel. If higher than predicted, it adds the value (number of observed counts)–(number of predicted counts) to the list; otherwise, it adds a zero. Finally, the sum of the values for each photon is normalized to 1, and the obtained list of probabilities is added to the input FITS file. We note that the "other = yes" mode can be used only to search for new sources, but not for the analysis of an already-detected source: in fact, this would induce a systematic error due to the wrong characterization of the statistics. On the other hand, accurate analysis of the residuals would allow one to find low-statistic point-like sources or extended sources that can be added to the input list to iteratively improve the results of the script. In order to extract photons from the source of interest, only sources ∼75 away from it should be considered. The XMM-Newton theoretical PSF is null farther than 5' from the source, so that only sources that affect this circle are necessary.

In order to validate this method, we applied it to several different pulsar fields: we used the weighted MCMC analysis presented in Section 4 to compare the H-value resulting from weighted and unweighted tests. We obtained improvements of H-value ranging from 10% to 50%, with better results for pulsars with multicomponent spectra, superimposed on nebular and/or point-like sources. In the case of the already-published J1813–1246 (Marelli et al. 2014b), a simpler version of this script increased the H-value from 11123 to 12092; the use of this new method further increased it to 12432. In the case of J2055, the H-value increased from 52 to 66 (with the same number of trials). In the case of Geminga, we used some of the public XMM-Newton observations to detect the nonthermal pulsation, in the 2–10 keV band. Using our new script, we found an increase in the H-value from 120 to 220, almost doubling it. Tests on other pulsars are in progress. We also considered the distribution of residuals of the considered point-like sources. For each of the ∼50 objects considered, we found residuals of no more than a few percent of the original fluxes, thus confirming the validity of the method.

The case of source 1 (see Table 1) is representative of the capacity of this script. That object was revealed by both the SAS tool edetect_chain and the CIAO tool wavdetect as a single point-like emission. A brightness profile analysis does not reveal significant deviations from the theoretical XMM-Newton PSF. A spectral analysis gives instead a peculiar result: this source cannot be fitted by a single-model spectrum but requires a double-component spectrum (apec+apec or apec+power-law). Our script allowed for a more accurate analysis of the source. Taking into account a single point-like object as source 1, we obtain residuals of more than 40% of source counts (while the other ∼50 sources give only a few percentage points), distributed in pixels southeast to northwest from the best source 1 position; we also obtain an overestimation of counts in pixels southwest to northeast, so that the total number of simulated counts is in agreement with the observed ones. We note that the nearby nebula should not significantly contribute to the source counts, owing to its faintness. Then, we tried a model with two different sources: the positions were iteratively refined in order to minimize the residuals. We found that two sources separated by ∼8'' describe adequately the "source 1 system," with lower residuals of ∼15%. The brightest one has a flux ∼10 times larger than the weakest, and its spectrum is best described by a power law, while the other one spectrum is best described by an apec; spectral parameters and optical counterparts allow us to associate them with an AGN and a star, respectively. We note that a three-source system, forming an equilateral triangle with a side of ∼10'', would lower the residuals to ∼3%.

Although this program is much more evolved than the script we used for the analysis of PSR J1813–1246 (Marelli et al. 2014b), it still needs further technical improvement to reduce the high computing time. A future version of this script could also be more user-friendly: a simplified version of the analysis could be run using the energy bands and associated CRs for each of the objects reported in the 3XMM catalog of serendipitous XMM-Newton sources.¹¹ Future implementations could also allow the user to analyze data from X-ray observatories different than XMM-Newton, e.g., NuSTAR (Harrison et al. 2013), by simply changing the PSF model and some variables.