A DEEP CAMPAIGN TO CHARACTERIZE THE SYNCHRONOUS RADIO/X-RAY MODE SWITCHING OF PSR B0943+10

PSR B0943+10 was observed with the LOFAR (Stappers et al. 2011; van Haarlem et al. 2013) for a total of ∼56 hr. The LOFAR international stations in the United Kingdom, France, Sweden, and Germany (stations UK608, FR606, SE607, and DE601–DE605) were used because these provide adequate sensitivity for separating PSR B0943+10's B- and Q-modes and are less oversubscribed than the central part of the LOFAR interferometric array, which is located in the Netherlands. Also, having multiple stations record simultaneously greatly lowers the risk that there will be a gap in the radio coverage during the XMM-Newton observations.

Each LOFAR international station includes 96 highband antenna (HBA) tiles, which are sensitive to radio waves in the range 110–250 MHz (though not the whole frequency range at once). The signals from individual HBA elements are coherently summed at the station, producing a combined field of view that is effectively like that of a single-dish radio telescope. The station signals from UK608 were recorded locally with an Artemis system (Karastergiou et al. 2015). The signals from SE607 and FR606 were recorded with the LuMP software²² on local recording machines, while the signals from the German stations were transmitted via the German LOFAR fiber network to the MPIfR in Bonn and recorded with LuMP there. For each station, we recorded 95.3 MHz of bandwidth (using all 488 "beamlets," of 195.3125 kHz each, in the LOFAR 8 bit station mode), at a center frequency of 149.9 MHz.

Initial data analysis and excision of radio-frequency interference (RFI) was done using the PSRCHIVE suite of pulsar software (van Straten et al. 2012), and associated scripts from Coastguard ²³ (Lazarus et al. 2016). We folded these data sets using an up-to-date rotational ephemeris derived from ongoing monitoring of PSR B0943+10 using the Lovell telescope (see Table 4), but using a contemporaneous dispersion measure derived from the LWA data sets to de-disperse the data.

The Arecibo observations were conducted using the 305 m telescope and 327 MHz receiver system in Puerto Rico. The pulse sequences were acquired using four Mock spectrometers covering adjacent 12.5 MHz portions of the total 50 MHz band. Arecibo rise-to-set observations of about 2.4 hr or, equivalently, 8000 pulses were carried out for each session apart from that of November 21, which was truncated by scheduling exigencies. These observations were processed into calibrated polarimetric pulse sequences with milli-period time resolution, wherein the modal transitions could readily be identified by visual inspection and then confirmed by fluctuation-spectral analysis (apart from the final session wherein persistent RFI made identification difficult).

The LWA data were taken using the first station of the LWA (LWA1: Taylor et al. 2012; Ellingson et al. 2013), which consists of 256 dipole antennas sensitive to radio waves in the range 10–88 MHz. LWA1 is co-located with the Very Large Array in New Mexico, USA. PSR B0943+10 was observed using LWA1 for a total of ∼72 hr in raw beam-forming mode using two beams (each having two tunings) in "split bandwidth" mode with center frequencies of 57.15, 64.5, 71.85, and 79.2 MHz. The observation setup was similar to other LWA1 pulsar observations described in Stovall et al. (2015). However, in order to reduce the data size due to the long observing duration, data were recorded with a reduced sample rate of 9.8M samples ${{\rm{s}}}^{-1}$. The resulting data were then coherently de-dispersed, divided into 256 frequency channels, and folded into 1024 profile bins with 60 s integrations using the dspsr software package (van Straten & Bailes 2011). RFI was excised using a median zapping routine from PSRCHIVE (van Straten et al. 2012) followed by manual zapping of RFI. The data were reduced to 16 frequency channels for subsequent analysis.

PSR B0943+10 was observed in the first half of each XMM-Newton observation by LOFAR and in the second half by LWA. Arecibo observed near the midpoint of the sessions and provided an important bridge in sensitivity when PSR B0943+10 was low on the horizon for both LOFAR and LWA (see Table 1). We combined the radio data into pulse profile stacks as a function of time, using the PSRCHIVE archive format to store and manipulate the data. The individual pulses were summed in 5 minute blocks of time to increase the signal-to-noise ratio (S/N). Figure 1 shows a representative sample of the data for Session 5 alongside the corresponding XMM-Newton X-ray light curve. The radio profile stacks for all observing sessions are shown in the Appendix. Each 5 minute integration was qualified as being either B- or Q-mode, based on the S/N and pulse morphology. We note that the pulse profile morphology of PSR B0943+10 evolves rapidly at low radio frequencies (Bilous et al. 2014). Therefore, the LWA (40–80 MHz), LOFAR (110–190 MHz), and Arecibo (302–352 MHz) profile morphologies are different. Table 2 provides a summary of all modal time ranges and their durations. The shortest observed mode instance lasted 22 minutes (B-mode at the beginning of the November 23 observation, but this is a lower limit to the actual duration since the previous mode transition was not observed). The longest instance was longer than 17 hr (Q-mode of the last observation); this is, to our knowledge, the longest mode duration ever reported for this pulsar. We have determined the mode instances using a time resolution of 5 minutes because the S/N of the profiles using shorter integrations is arguably too low to make a robust mode identification. We note that while some mode-changing pulsars have in some cases shown mixing between the modes, making any such identification dubious, this is not the case for PSR B0943+10. Given the low X-ray count rate of PSR B0943+10 (at most ∼0.011 counts s⁻¹), and the fact that the modes typically last hours, the uncertainty in the times of the mode transitions due to the adopted 5 minute resolution is negligible.

Zoom In Zoom Out Reset image size

The pn and MOS counts for the timing analysis were extracted from a circular region with a radius of 15'', and the times of arrival were converted to the barycenter of the solar system. We excluded the parts of the observations affected by a high background level by selecting only time intervals in which the EPIC pn count rate in the range 10–12 keV was below 1.2 counts s⁻¹. This resulted in exposure times of 119 ks and 175 ks for the Q- and B-modes, respectively.

The pulse phases of the PSR B0943+10 counts were computed using the same ephemeris adopted for the radio data (Table 4). In Figure 2 we show the folded pulse profiles in different energy ranges for the Q- and B-mode time intervals, as well as those corresponding to the whole observation. The latter differs from the sum of the B and Q time intervals, since it includes also ∼17 ks for which information on the radio mode is missing.

Zoom In Zoom Out Reset image size

Table 4.  Ephemeris of PSR B0943+10 Derived from Monitoring Observations with the Lovell Telescope

R.A. (J2000)	09h46m07786
Decl. (J2000)	09° 52' 0076
ν (Hz)	0.9109890963(2)
$\dot{\nu }$ (Hz s−1)	−2.9470(7) × 10−15
$\ddot{\nu }$ (Hz s−2)	−1.38 × 10−25
Epoch (MJD)	55962.057374
Validity range (MJD)	54861.014–57011.249

Note. The coordinates refer to the X-ray position of PSR B0943+10, which has an accuracy of $1^{\prime\prime}$ (Hermsen et al. 2013). Numbers in parenthesis are the 1σ errors on the least significant digits.

Download table as:  ASCIITypeset image

The pulsations are clearly detected in the energy range 0.5–2 keV during both radio modes, with a statistical significance, computed using the Rayleigh test statistics ${Z}_{1}^{2}$ (Buccheri et al. 1983), of 8.4σ in the Q-mode and 6.3σ in the B-mode. The 0.5–2 keV pulse profiles of the two modes are similar. They are broadly sinusoidal, with pulsed fractions, defined as the amplitude of the sinusoid divided by the average value, of (45 ± 6)% (Q-mode) and (47 ± 6)% (B-mode), and aligned in phase. These pulsed fractions are consistent with the more precise values obtained in Section 3.3 with a maximum likelihood (ML) analysis. By a cross-correlation of the two light curves, we could set an upper limit at a 90% confidence level of 0.08 cycles on the relative phase of the peak in the B- and Q-modes. No pulsations are seen above 2 keV, while there is a hint that they could be present also below 0.5 keV. In fact, although the statistical significance of the modulation in the range 0.2–0.5 keV is only $\sim 2\sigma$ in the Q- and B-modes separately, a level of 3.4σ is reached when one considers the whole data set.

As a first step in the analysis of the phase-averaged spectra of PSR B0943+10, we reprocessed the pn data using the SAS task epreject to reduce the detector noise at the lowest energies. We then applied a tight screening to eliminate periods of high background due to flares of soft protons. For the extraction of the source spectra we used circular regions of radius 15'', while for the background we used circles of radius 35'' for the pn and 25'' for the MOS. We used single-pixel events for the pn (PATTERN = 0) and single- and multiple-pixel events for the MOS (PATTERN ≤ 12). We first combined the spectra of the two MOS cameras and then summed the spectra of the Q- and B-mode time intervals of all the observations. This resulted in four spectra with net live times of 67.7 ks (pn) and 113.5 ks (MOS) in the Q-mode, and of 90.9 ks (pn) and 162.6 ks (MOS) in the B-mode.

For each mode, we simultaneously fitted the pn and MOS spectra in the energy range 0.2–10 keV using XSPEC. Following previous works (Zhang et al. 2005; Hermsen et al. 2013) and to facilitate comparisons with their results, we fixed the absorbing column density to ${N}_{{\rm{H}}}=4.3\times {10}^{20}$ cm⁻². This value corresponds to the dipersion measure of PSR B0943+10 assuming a 10% ionization of the interstellar medium (e.g., He et al. 2013).

The analysis of the Q-mode spectra showed that fits with a single power law or a single blackbody are not acceptable. With fixed N_H, these models gave values of ${\chi }_{\nu }^{2}\gt 2.5$ (for 20 degrees of freedom (dof)), corresponding to null hypothesis probabilities (nhp) smaller than 2 × 10⁻⁴. Letting N_H vary freely, the single-blackbody fit was still rejected (nhp = 2 × 10⁻³), while a formally acceptable fit was obtained with a power law (${\chi }_{\nu }^{2}$/dof = 1.4/19, nhp = 0.12). However, the latter resulted in a large power-law photon index (${\rm{\Gamma }}={3.3}_{-0.2}^{+0.3}$) and an absorption value (${N}_{{\rm{H}}}=(1.8\pm 0.4)\times {10}^{21}$ cm⁻²) much higher than that expected for PSR B0943+10 (the total Galactic column density in this direction is $2.5\times {10}^{20}$ cm⁻², Kalberla et al. 2005). A good fit to the Q-mode spectra was instead obtained by using either a blackbody plus power-law model or the sum of two blackbodies, resulting in the parameters given in Table 5.

Table 5.  Results for the Total Spectra of Q- and B-modes with Standard Analysis

	Q-mode	Q-mode	B-mode	B-mode	Joint fit of Q- and B-modes
	BB + PL	BB + BB	BB	BB + PL	BB + PL
Photon index	${2.7}_{-0.2}^{+0.3}$	⋯	⋯	2.5 ± 0.4	${2.6}_{-0.2}^{+0.3}$
${K}_{\mathrm{PL}}$ a	2.4 ± 0.6	⋯	⋯	0.6 ± 0.5	2.1 ± 0.5b/0.9 ± 0.2c
${F}_{\mathrm{PL}}$ d	5.5 ± 1.1	⋯	⋯	1.4 ± 0.4	⋯
kT1 (keV)	0.29 ± 0.03	0.14 ± 0.02	0.23 ± 0.01	0.24 ± 0.03	0.27 ± 0.02
${\alpha }_{\mathrm{BB}1}$ e	${1.15}_{-0.40}^{+0.75}$	${29}_{-10}^{+18}$	${2.2}_{-0.4}^{+0.6}$	${1.4}_{-0.5}^{+1.1}$	1.6 ± 0.6b/0.7 ± 0.4c
${R}_{\mathrm{BB}1}$ f (m)	${21}_{-4}^{+6}$	${106}_{-20}^{+29}$	${29}_{-3}^{+4}$	${23}_{-5}^{+8}$	25 ± 5b/16 ± 4c
${F}_{\mathrm{BB}1}$ d	6.4 ± 1.2	4.8 ± 1.3	⋯	${3.8}_{-0.8}^{+2.4}$	⋯
kT2 (keV)	⋯	0.37 ± 0.04	⋯	⋯	⋯
${\alpha }_{\mathrm{BB}2}$ e	⋯	${0.51}_{-0.20}^{+0.33}$	⋯	⋯	⋯
${R}_{\mathrm{BB}2}$ f (m)	⋯	${14}_{-3}^{+4}$	⋯	⋯	⋯
${F}_{\mathrm{BB}2}$ d	⋯	7.5 ± 1.3	⋯	⋯	⋯
${F}_{\mathrm{TOT}}$ d	11.8 ± 0.5	12.2 ± 0.5	5.2 ± 0.3	5.2 ± 0.3	⋯
${\chi }_{\nu }^{2}$/dof	1.24/18	0.98/18	1.12/20	1.03/18	1.11/39
nhp	0.22	0.47	0.32	0.42	0.29

Notes. Joint fits of pn + MOS spectra with fixed ${N}_{{\rm{H}}}=4.3\times {10}^{20}$ cm⁻² (BB = blackbody, PL = power law). Errors at 1σ.

^aNormalization of the power law at 1 keV in units of 10⁻⁶ photons cm⁻² s⁻¹ keV⁻¹. ^bQ-mode. ^cB-mode. ^dFlux in the range 0.5–2 keV corrected for the absorption, in units of 10⁻¹⁵ erg cm⁻² s⁻¹. ^eBlackbody normalization in units of 10⁻⁴ photons cm⁻² s⁻¹ keV⁻³. ^fBlackbody radius for d = 630 pc.

Download table as:  ASCIITypeset image

The B-mode spectra were well fit by a blackbody with temperature kT = 0.23 keV, while a power law was clearly rejected (${\chi }_{\nu }^{2}$/dof = 2.1/20, nhp = 3 × 10⁻³). A single power law gave a statistically acceptable fit (nhp = 0.13) but with rather implausible parameters (${\rm{\Gamma }}={4.1}_{-0.6}^{+0.8}$, ${N}_{{\rm{H}}}=(3.2\pm 0.1)\times {10}^{21}$ cm⁻²). Good fits could also be obtained with two-component models, but, contrary to the case of the Q-mode, the addition of a further component to the single blackbody is not statistically required.

Finally, to assess whether there is a significant spectral difference between the two modes, we fitted simultaneously the B and Q spectra with a blackbody plus power-law model keeping Γ, kT, and relative normalization of the two components tied to common values. In this way we obtained a good fit with the parameters given in the last column of Table 5 and a scale factor between the B and Q spectra of ${f}_{{\rm{B}}/{\rm{Q}}}=0.42\pm 0.03$.

In this subsection we follow an alternative approach, based on an ML analysis, which is particularly powerful in the case of faint sources. This method (2D-ML, in the following) has been used by Hermsen et al. (2013) to obtain the source count rates and spectra in Q- and B-modes for the first radio/X-ray campaign on PSR B0943+10 performed in 2011. Here we extend it to derive also the spectra of the pulsed and unpulsed components of both radio modes (3D-ML).

Contrary to the traditional analysis, in which the background is estimated from a "source-free" region of the image, in the 2D-ML the instrumental point-spread function (PSF) and a background distribution are simultaneously fitted to the data in the region including the source of interest, taking into account the Poissonian nature of the counting process. At first, the selected events are sorted according to their spatial coordinates (x, y) to produce skymaps of counts.²⁴ The likelihood function L is defined as

$$\begin{eqnarray*}L & = & \mathrm{ln}({\prod }_{i,j}\ ({\mu }_{{ij}}^{{N}_{{ij}}}\exp (-{\mu }_{{ij}}){/N}_{{ij}}!))\\ & = & \displaystyle \sum _{i,j}{N}_{{ij}}\mathrm{ln}({\mu }_{{ij}})-{\mu }_{{ij}}-\mathrm{ln}({N}_{{ij}}!),\end{eqnarray*}$$

where ${\mu }_{{ij}}=\beta +\sigma \times {\mathrm{PSF}}_{{ij}}$ is the expectation value for pixel (i, j) and N_ij is the number of counts measured in pixel (i, j). The quantity L is maximized simultaneously with the background parameter β and the source scale parameter σ. Since the PSF is normalized to unity, σ gives directly the total number of net (background-free) source counts. The second derivative matrix of L, evaluated at the maximum value, contains information on the uncertainties of the derived parameters, β and σ. Applying this optimization scheme for a grid of energy windows results in net source counts per energy interval, which can be converted to source fluxes by a forward-folding fitting procedure by assuming a spectral model and using appropriate response (arf and rmf files) and livetime information (dead-time-corrected exposure).

We applied this method, separately for the pn and for the sum of the two MOS CCDs, to extract the time-integrated Q- and B-mode spectra in the range 0.2–10 keV. Time intervals with high background were excluded using the same cuts of Section 3.1. For the pn and MOS PSF we used axially symmetric profiles as derived from in-flight calibrations²⁵ , with parameters appropriate for an energy of 1.5 keV. We modeled the background with a uniform spatial distribution. As in the previous subsection, we fixed ${N}_{{\rm{H}}}=4.3\times {10}^{20}$ cm⁻².

For the Q-mode, single-component power-law or blackbody models yielded unacceptable ${\chi }^{2}$ values, while for the B-mode a single power law is excluded, contrary to the 2011 data, which could be described equally well by either a power law or a blackbody. A single blackbody, however, adequately describes the data, as does the combination of a blackbody plus either a power law or another blackbody. These results, summarized in Table 6, are fully consistent with those described in the previous subsection.

Table 6.  Results for the Total Spectra of Q- and B-modes with 2D-ML Analysis

	Q-mode	Q-mode	B-mode	B-mode	Joint fit of Q- and B-modes
	BB + PL	BB + BB	BB	BB + PL	BB + PL
Photon index	2.4 ± 0.2	⋯	⋯	${2.3}_{-0.5}^{+0.4}$	2.4 ± 0.1
KPL	3.2 ± 0.7	⋯	⋯	0.8 ± 0.5	2.8 ± 0.2a/1.2 ± 0.2b
FPL	7.4 ± 0.7	⋯	⋯	1.6 ± 0.4	⋯
kT1 (keV)	0.27 ± 0.04	0.11 ± 0.02	0.24 ± 0.01	0.24 ± 0.03	0.26 ± 0.01
${\alpha }_{\mathrm{BB}1}$	${1.12}_{-0.49}^{+0.85}$	${61}_{-32}^{+77}$	${2.3}_{-0.4}^{+0.5}$	${1.5}_{-0.6}^{+1.0}$	1.7 ± 0.2a/0.8 ± 0.1b
RBB1 (m)	${21}_{-5}^{+8}$	${153}_{-40}^{+97}$	${29}_{-3}^{+3}$	${24}_{-5}^{+8}$	25 ± 2a/17 ± 1b
FBB1	4.9 ± 0.9	3.0 ± 0.4	⋯	4.1 ± 0.6	⋯
kT2 (keV)	⋯	0.34 ± 0.03	⋯	⋯	⋯
${\alpha }_{\mathrm{BB}2}$	⋯	${0.90}_{-0.26}^{+0.36}$	⋯	⋯	⋯
RBB2 (m)	⋯	${19}_{-3}^{+4}$	⋯	⋯	⋯
FBB2	⋯	9.6 ± 0.9	⋯	⋯	⋯
FTOT	12.3 ± 1.1	12.6 ± 1.0	5.7 ± 0.3	5.7 ± 0.7	⋯
${\chi }_{\nu }^{2}$/dof	0.63/30	1.02/30	1.15/24	0.85/22	0.71/56
nhp	0.94	0.44	0.28	0.66	0.95

Notes. Joint fits of pn + MOS spectra with fixed ${N}_{{\rm{H}}}=4.3\times {10}^{20}$ cm⁻². Errors at 1σ.

^aQ-mode. ^bB-mode.

Download table as:  ASCIITypeset image

Figure 4 shows with black symbols (circles for pn, squares for the summed MOS) the unabsorbed Q-mode (left) and B-mode (right) spectra for the best-fit blackbody plus power-law model. The separate spectral components are shown as black dashed lines, while the black solid lines are their combination.

As first done in Hermsen et al. (2016) for the mode-switching pulsar PSR B1822−09, we can easily generalize the 2D-ML method to take into account also the pulse phase information of the events. In this 3D-ML approach, we can sort the events according to their spatial coordinates and pulse phase (x, y, ϕ). The expectation value of bin (i, j, k) can now be written as

$$\begin{eqnarray*}&&{\mu }_{{ijk}}=\beta +{\sigma }_{u}\times {\mathrm{PSF}}_{{ij}}+{\sigma }_{p}\times {\mathrm{PSF}}_{{ij}}\times {{\rm{\Phi }}}_{k}.\end{eqnarray*}$$

Here the value of the normalized pulse profile at bin k is represented by ${{\rm{\Phi }}}_{k}$, while ${\sigma }_{u}$ and ${\sigma }_{p}$ correspond to the scale factors for the unpulsed and pulsed components, respectively. The pulsed fraction η can be determined as $\eta =1/[1+({N}_{\phi }\times {\sigma }_{u}/{\sigma }_{p})]$, with ${N}_{\phi }$ the number of bins of the normalized pulse profile. As shown in Section 3.1, the X-ray pulse profile of PSR B0943+10 is well described by a sinusoid with the same phase in both radio modes. We assume in the following an energy-independent sinusoidal pulse profile.

With this method we derived the pulsed fraction as a function of energy for the Q- and B-modes shown in Figure 3. The Q-mode pulsed fraction increases steadily from ∼0.21 in the 0.2–0.5 keV band to ∼0.69 at 2 keV (consistent with the values reported by Hermsen et al. 2013), and then drops significantly. This behavior is different from that shown in the B-mode, in which the pulsed fraction is constant across the 0.2–2 keV energy band, with an average 0.5–2 keV value of (38 ± 5)% and even a significant measurement of (35 ± 11)% in the 0.2–0.5 keV band. The 0.5–2 keV pulsed fraction for the Q-mode is (52 ± 4)%.

Zoom In Zoom Out Reset image size

We also performed a spectral analysis of the unpulsed and pulsed components of both modes using the source counts per energy bin extracted with the 3D-ML approach. The results are summarized in Table 7. For the Q-mode (10 spectral bins in the 0.2–10 keV band) we found that the pulsed component is properly described by a blackbody (Figure 4(a); red data points; red solid line for model), while a power law is rejected (${\chi }_{\nu }^{2}/\mathrm{dof}=5.09/8$). On the other hand, the unpulsed component (Figure 4(a); blue data points; blue solid line for model) is well fit by a power law, while a blackbody is rejected (${\chi }_{\nu }^{2}/\mathrm{dof}=3.53/8$). The sum of the two fits is indicated by a purple solid line in Figure 4(a), which coincides within uncertainties with the total spectrum derived with the 2D-ML method (black line).

Zoom In Zoom Out Reset image size

Table 7.  Results for the Pulsed and Unpulsed Spectra of Q- and B-modes with 3D-ML Analysis

	Q-mode	Q-mode	B-mode	B-mode	B-mode	B-mode
	Pulsed	Unpulsed	Pulsed	Pulsed	Unpulsed	Unpulsed
	BB	PL	BB	PL	BB	PL
Photon index	⋯	2.5 ± 0.2	⋯	2.3 ± 0.3	⋯	2.4 ± 0.2
KPL	⋯	2.9 ± 0.3	⋯	0.9 ± 0.1	⋯	1.4 ± 0.2
kT (keV)	0.29 ± 0.02	⋯	${0.22}_{-0.03}^{+0.04}$	⋯	${0.22}_{-0.02}^{+0.03}$	⋯
${\alpha }_{\mathrm{BB}}$	${1.11}_{-0.28}^{+0.39}$	⋯	${1.24}_{-0.64}^{+1.24}$	⋯	${1.88}_{-0.72}^{+1.13}$	⋯
RBB (m)	20.6 ± 3.1	⋯	21.8 ± 8.3	⋯	26.9 ± 6.6	⋯
F	${6.1}_{-0.5}^{+0.6}$	${6.5}_{-0.5}^{+0.6}$	2.1${}_{-0.3}^{+0.4}$	2.0 ± 0.3	3.6 ± 0.4	3.2 ± 0.4
${\chi }_{\nu }^{2}$/dof	1.07/8	0.85/8	0.57/5	0.37/5	0.56/5	0.86/5
nhp	0.38	0.56	0.72	0.87	0.73	0.51

Note. Joint fits of pn + MOS spectra with fixed ${N}_{{\rm{H}}}=4.3\times {10}^{20}$ cm⁻². Errors at 1σ.

Download table as:  ASCIITypeset image

For the B-mode the characterization of the spectra of the unpulsed and pulsed components is more difficult due to the lower number of events (by ∼35%) of the pulsed and unpulsed emission. Therefore, we used only seven spectral bins. The unpulsed and pulsed components can be fit equally well by either a blackbody or a power law (see the last four columns of Table 7). Figure 4(b) shows the spectra of the pulsed and unpulsed emissions of the B-mode, assuming the same characteristics found for the Q-mode, i.e., thermal pulsed emission and non-thermal unpulsed emission; the red data points and solid line represent the (unabsorbed) pulsed flux measurements and best-fit blackbody model, respectively, and the blue data points and solid line the unpulsed flux measurements and best-fit power-law model. The purple solid line is the sum of the two models. The match with the 2D-ML results is less obvious than for the case of the Q-mode (note that the error regions are substantial), which is not strange given the degenerate nature of the B-mode spectrum.

Radio studies showed that, while the Q-mode emission is steady and largely chaotic, the highly organized B-mode evolves in its profile and subpulse behavior (Rankin & Suleymanova 2006). The evolution of the radio properties begins immediately after the Q to B transition and is largely accomplished within two hours, as can be seen in the radio charts displayed in the Appendix. Motivated by these results, we searched for a possible evolution of the X-ray properties.

We divided all the B-mode data into two subsets: one with all the data collected within 3 hr after the start of the B-mode (early-B in the following) and one with the remaining data (late-B). We excluded the initial 3 hr of the observations that started with the pulsar in B-mode but for which the transition time was not known. Although the radio evolution occurs on a shorter timescale (Bilous et al. 2014), our choice of 3 hr was dictated by the need for comparable counting statistics in the two subsets (351 and 362 pn+MOS counts in the early-B and late-B intervals, respectively).

Analyzing the data as in Section 3.1, we found that the pulsations in the range 0.5–2 keV are detected with a significance of only 2.6σ in the early-B subset, while their significance in the late-B subset is 4.6σ. By means of Monte Carlo simulations, assuming the average pulsed fraction of the whole B-mode data derived in Section 3.1, we found that the probability of obtaining a significance as small as that observed in the early-B data is ∼6%. The pulse profiles of the two subsets, plotted in Figure 5, show a different degree of modulation, also confirmed by the pulsed fractions obtained with the 3D-ML analysis: (27 ± 8)% for the early-B and (42 ± 8)% for the late-B. These results suggest that there might be an evolution of the X-ray properties during the B-mode, but the statistics are too low to draw firm conclusions. A similar analysis to search for variations in the pulsed fraction during the Q-mode gave negative results.

Zoom In Zoom Out Reset image size

To search for possible evolution in the X-ray flux, we also constructed stacked light curves of the mode transitions by summing all the data. These are plotted in Figure 6, where we have set the origin of the time axis to the time of the mode transitions. The count rates before and after the transitions are well fit by constant functions and there is no evidence for any gradual increase/decrease of the flux leading to, or following, the mode-switch times. The changes in the X-ray flux occur on a timescale shorter than the bin time of 900 s used in these light curves.

Zoom In Zoom Out Reset image size

To search for the presence of diffuse emission around PSR B0943+10 on angular scales smaller than those resolved with XMM-Newton, we carried out an observation with the Chandra satellite on 2016 January 15 (Obs. ID 16759). The observation was done with the ACIS-S instrument in Timed Exposure full-imaging mode (frametime: 3.14 s) and lasted 50 ks. The data reduction and analysis were perfomed with the CIAO analysis software v. 4.8 and the PSF simulation package ChaRT/MARX (v. 5.2), using the calibration files in the CALDB database v. 4.7.1.

PSR B0943+10 was clearly detected with a net count rate of $(1.0\pm 0.1)\times {10}^{-3}$ counts s⁻¹ in the energy range 0.3–8 keV. We extracted a spectrum using the counts in a circle of radius 125 (encircled energy fraction: EEF ∼ 95%) for the source and an annulus with radii 5'' and 10'' for the background. The spectrum is well fit by a blackbody model with absorption and temperature fixed to the values derived with XMM-Newton for the B-mode spectrum. The resulting flux is $(4.7\pm 0.8)\times {10}^{-15}$ erg cm⁻² s⁻¹ (0.5–2 keV), consistent with that of the B-mode. A fit with the blackbody plus power-law model with parameters fixed to those of the Q-mode is unacceptable, unless a rescaling by a factor 0.47 ± 0.07 is applied. Therefore, we conclude that PSR B0943+10 was in the B-mode during most (or possibly the whole) of the Chandra observation.

The radial distribution of the source counts around the pulsar position is consistent with the instrumental PSF, with no evidence for diffuse emission at radii larger than ∼2'' (EEF approaching 100%). We derived 3σ upper limits of 3 × 10⁻⁶ and $1.5\times {10}^{-6}$ counts s⁻¹ arcsec⁻² on the average surface brightness within radii of 5'' and 10'', respectively. For a power-law spectrum with photon index 2.5, these values correspond to flux upper limits of $1.0\times {10}^{-15}$ and $2.3\times {10}^{-15}$ erg cm⁻² s⁻¹ (0.5–2 keV). The limits derived for larger nebulae are less constraining: for example, in the case of diffuse emission extending to 20'', the upper limit on surface brightness is 10⁻⁶ counts s⁻¹ arcsec⁻², corresponding to 7 × 10⁻¹⁵ erg cm⁻² s⁻¹.

The simultaneous emission of thermal and non-thermal X-rays is a common property seen in all the rotation-powered pulsars that have been studied with sufficient sensitivity. The thermal emission from young and middle-aged neutron stars results from internal cooling and involves a relatively large fraction of (or the whole) surface of the star. In older objects, such as PSR B0943+10, only small regions at the magnetic poles are sufficiently heated by backward-accelerated magnetospheric particles to significantly emit in the X-ray band.

The presence of thermal and non-thermal X-rays during the Q-mode was already clear in the 2011 data and is confirmed by the new observations, which yield fully compatible spectral parameters. Furthermore, the analysis of the Q-mode data with the powerful 3D-ML method shows that the pulsed component is well fit by a blackbody and not by a power law, as found by Hermsen et al. (2013). The opposite is true for the unpulsed component, which is clearly inconsistent with a blackbody. Therefore, we can describe the Q-mode emission as the sum of two contributions with similar fluxes in the range 0.5–2 keV: a pulsed thermal component, with blackbody temperature kT ∼ 0.28 keV and emission radius ∼21 m, plus an unpulsed non-thermal component, with power-law photon index ∼2.5.

Although a single blackbody, with kT = 0.23 keV and emission radius ∼30 m, gives a fully satisfactory fit of the B-mode spectrum, the presence of an additional power law is fully compatible with the data. This yields a smaller emitting area for the blackbody, similar to that seen in the Q-mode. The 3D-ML analysis indicates that in the B-mode both the pulsed and unpulsed spectra can be fit equally well by either a power law or a blackbody (see Table 7). However, from a physical point of view, it is quite contrived to invoke a mechanism able to rapidly switch off (or strongly suppress) the pulsed thermal component seen in the Q-mode and replace it with a non-thermal pulsed emission. It is much more natural to explain the B-mode pulsed flux as thermal emission closely related to that seen in the Q-mode, and possibly of the same origin. This is also supported by the fact that the shape of the pulse profile does not change, being nearly sinusoidal with the maximum very close in phase to the main radio pulse in both modes.

It is interesting to note that the best-fit value of the temperature is slightly higher during the Q-mode. The 30% change in temperature corresponds, for a fixed emitting area, to an increase in bolometric luminosity by a factor ∼3, consistent with the measured flux increase of 2.9. This increase in the thermal flux can be compared to the factor ∼2 increase in the non-thermal luminosity. The error regions for the parameters of the blackbody components (Figure 7) are consistent with either a change in emission area (with fixed temperature) or a change in temperature (with fixed area), although the best-fit values favor the latter possibility.

Zoom In Zoom Out Reset image size

Drifting radio subpulses during the B-mode of PSR B0943+10 serve as a classical example for associating the phenomenon with the ${\boldsymbol{E}}$ × ${\boldsymbol{B}}$ drift of a system of sparks that are generated at the inner accelerating region (IAR). The theory was first suggested by Ruderman & Sutherland (1975), who conjectured the IAR to be an inner vacuum gap (IVG). If ${\boldsymbol{\Omega }}\cdot {{\boldsymbol{B}}}_{{\rm{s}}}$ < 0 (where ${\rm{\Omega }}=2\pi /P$, P being the pulsar period and B_s the surface magnetic field), an IVG can be formed just above the polar cap, where a high potential difference exists. Such a region can discharge as a system of sparks by the process of magnetic pair creation. The electric field in the IVG separates positrons and electrons. The latter accelerate back to the polar cap and heat the surface sufficiently enough to generate soft X-ray thermal emission. The positrons instead accelerate away from the stellar surface and produce secondary pairs, which eventually leads to non-stationary spark-associated plasma columns, where coherent radio emission is generated at a distance of about 50 stellar radii. The plasma columns undergo an ${\boldsymbol{E}}$ × ${\boldsymbol{B}}$ drift in the IAR, which is reflected in the radio emission region.

Deshpande & Rankin (2001) were able to model the drifting radio features observed in the B-mode of PSR B0943+10 as a system of 20 sparks rotating in the outer edge of the pulsar beam as a carousel, with a circulation time of 37P. The fact that this circulation time was longer than that predicted by Ruderman & Sutherland (1975) led to the refinement of the IAR model as a partially screened gap (PSG) (Gil et al. 2003). The basic feature of the PSG model is that the polar cap is maintained just above a critical temperature T_c where positive ions can be extracted from the surface, which can then screen the vacuum electric field and reduce the drift speed.

Our blackbody fit to the pulsed B-mode spectrum gives an emitting area A_s of ∼1500 m², two orders of magnitude smaller than that expected for the polar cap in a dipole geometry ($A\sim \pi {R}_{\mathrm{pc}}^{2}=\tfrac{2{\pi }^{2}{R}^{3}}{{cP}}=1.3\times {10}^{5}$ m², for a stellar radius R = 13 km). This suggests that only a small fraction of the polar cap is heated or that the magnetic field is non-dipolar at the surface of the polar cap. Indeed such non-dipolar fields are expected in pulsars, as has been shown by Gil & Mitra (2001) and Gil et al. (2002a, 2002b). If the field near the surface is not dipolar, the ratio of the surface field to dipolar field would be ${B}_{{\rm{s}}}/{B}_{{\rm{d}}}\sim 100$, by flux conservation, implying ${B}_{{\rm{s}}}\sim 4\times {10}^{14}$ G. Geppert et al. (2013) showed that Hall drift of the crustal magnetic field can generate such strong fields.

No clear drifting radio features are seen during the Q-mode, so no further analysis regarding the radio and X-ray connection can be established, although, once again, the X-ray data appear to strongly suggest a very small blackbody emitting area, similar to that of the B-mode.

The next important phenomenon is the change in the pulsar radio emission from B- to Q-mode, where the tip-over occurs over perhaps hardly more than a single period. In the context of inner gap models, Zhang et al. (1997) proposed that the mode switching is caused by transitions between two states in which inverse Compton scattering plays a major role in causing the inner gap breakdown. At the mode-switch, the gap height and Lorentz factor of the primary particles change rapidly. The transitions can occur when the pulsar surface temperature is close to a critical value T_c, which, for PSR B0943+10 and in the case of a multipolar magnetic field, is predicted by these authors to be T_c = 10⁶ K. Considering the uncertainties involved on both the theoretical and observational sides, this value is in good agreement with our results.

More recently, Szary et al. (2015) suggested that the effect could arise due to two states of the PSG, termed the PSG-off and PSG-on states. In the PSG-off state the IAR is initially vacuum, and pair cascades due to curvature radiation can operate. The backflowing particles then hit the surface and raise the temperature to T_c, where the ions are extracted. In the PSG-on state the surface is already at a temperature slightly above T_c and the dominant process of pair cascade is due to inverse Compton scattering. The ions that are now extracted from the star can completely screen the electric field and no further heating can take place. Thus the prediction is that in both PSG-off and PSG-on states the temperature of the polar cap is $\sim {T}_{{\rm{c}}}$. One also does not expect a change in the surface area of the polar cap between the PSG-on/off states. The change in the intensity of the X-ray blackbody appears to be rather perplexing in this regard, thus ruling out a simple application of the PSG-on/off model.

The presence of the Q-mode precursor in the radio might be invoked to find clues about some additional source of change in area or temperature needed to account for the higher blackbody flux. If the precursor emission is associated with the closed field lines, that could lead to an increased area, while if it is associated with an increased return current, that causes an increase in temperature.

It has been argued that mode changes can be seen as changes in state of the global magnetospheric configuration. Models of the magnetosphere can be built that are physically self-consistent solutions of the highly nonlinear equations (Goodwin et al. 2004; Timokhin 2006), including a range of configurations for the last closed field line and hence possible polar-cap boundaries. However, these are steady-state models and give no hint as to what triggers a change from one state to another. Nor do they give the timescale of the mode durations or—especially—explain why certain states are preferred over others (and why, it seems, only two?).

Radio observations of other pulsars give clues that mode-changing involves not only the polar cap but the whole magnetosphere. Notably, Kramer et al. (2006), Lyne et al. (2010), Camilo et al. (2012), and Lorimer et al. (2012) reported evidence for a relation between radio mode switching and changes in period derivative, suggesting abrupt changes in the rate of angular momentum loss caused by sudden changes in magnetospheric currents along the open field lines. Such changes in magnetospheric currents could presumably be responsible for switches in non-thermal magnetospheric emission as well as in emitted thermal X-rays due to enhanced heating or changes in hot-spot area. Furthermore, PSR B1822–09, a nearly orthogonal rotator in which we can observe both poles, has two distinct modes: one with a regular pulse modulation and one largely chaotic in a manner resembling PSR B0943+10 (Gil et al. 1994; Latham et al. 2012). The modulations are in phase at both poles and the modes change rapidly and nearly simultaneously at both poles on a timescale of minutes. This suggests interpole communication. Similar phase-locked pulse modulation between poles has also been found in PSR B1702–19 and PSR B1055–52 (Weltevrede et al. 2007, 2012).

It may be that communication between poles is possible only when a pulsar is highly inclined, but it may simply be that this effect can be observed only when both poles are visible. If communication is possible for any modulating or mode-changing pulsar, we may speculate that it is also occurring in PSR B0943+10. The obvious line of communication would be via the magnetic field lines at or close to the polar-cap boundary. If modulation with the observed precision of the B-mode were occurring at both poles it would challenge the view that the drifting subpulses are the product of ${\boldsymbol{E}}$ × ${\boldsymbol{B}}$ drift in a local multipole configuration (it being already hard to understand how sparks could form a perfect carousel in such an environment) and would require the presence of magnetosphere-wide drift similar to auroral phenomena (e.g., Wright 2003).

If the ${\boldsymbol{E}}$ × ${\boldsymbol{B}}$ drift is a global phenomenon, one might expect a precise short $\approx 2P$ modulation, such as that found in PSR B0943+10, only to be sustainable for any length of time in the near-symmetry of a nearly aligned pulsar magnetosphere. Thus, although the B-mode can be seen as a steady magnetospheric eigenstate, the specific field lines of communication between the poles would themselves gradually evolve as a result of the slight misalignment of the rotation and magnetic axes. Eventually the delicately maintained eigenstate would suddenly collapse and a chaotic Q-mode ensue until the B-mode eigenstate could assert itself again. In this picture the observed temperature difference between the modes might be a consequence of the mode change and its evolution, rather than its cause, and the same could be argued for the non-thermal emission.

The 3D-ML analysis has shown that the pulsed component in the Q-mode is thermal, and the same conclusion may also hold for the B-mode, although with less certainty. The high pulsed fraction of this component is rather puzzling, since the radio emission properties have long suggested that PSR B0943+10 is an aligned rotator, and hence the polar-cap hot spot should not give rise to strongly pulsed X-rays.

To address this issue in some quantitative detail we computed the radiation intensity $I(E,\mu )$ emerging from a plane-parallel, magnetized atmosphere with ${\boldsymbol{B}}$ along the vertical axis, as is appropriate if the hot cap is at the magnetic pole; here E is the photon energy at the surface of the star and μ is the cosine of the angle between the surface normal and the ray direction. The magnetic field strength is $B=4\times {10}^{12}\,{\rm{G}}$ and we take an effective temperature, at infinity, of ${T}_{\mathrm{eff}}=250\,\mathrm{eV};$ the star's mass and radius are $M=1.4\ {M}_{\odot }$ and $R=13\,\mathrm{km}$, respectively. The composition of the atmosphere is pure H and complete ionization is assumed (see Zane & Turolla 2006, and references therein). $I(E,\mu )$ is then used to derive the pulse profiles, by means of the technique described in Turolla & Nobili (2013), which accounts for general-relativistic ray-bending and applies to circular spots of arbitrary size. Actually, since here the cap radius is $\ll R$, the emitting spot is virtually point-like and the use of a single atmospheric model is fully justified. In order to make the comparison with observational data meaningful, in computing the pulsed fraction the contribution of an unpulsed power-law component with the appropriate spectral parameters was taken into account.

The light blue region in Figure 3 shows the computed values of the pulsed fraction for the most likely ranges of the two geometrical angles, $5^\circ \lt \chi \lt 10^\circ$ and $10^\circ \lt \xi \lt 15^\circ$ (Deshpande & Rankin 2001); here ξ is the angle between the magnetic and rotation axes, and χ is the angle between the rotation axis and the line of sight. The computed pulsed fraction increases with energy then drops above ∼2 keV where the power-law flux becomes dominant, but it is systematically too low to account for the observed values, especially in the Q-mode. In a similar computation, Storch et al. (2014) adopted a model atmosphere with ${T}_{\mathrm{eff}}\sim 150\,\mathrm{eV}$. We recomputed the pulsed fraction following the same approach outlined above, but with an effective temperature of $140\,\mathrm{eV}$. While the pulsed fractions in this case are somewhat higher (darker region in Figure 3), they still fail to account for the observed ones.

The pulsed fractions computed by Storch et al., corrected for the presence of an unpulsed component, are closer to the data, although they still underpredict the observed values (see the dashed line in Figure 3). The main difference between the two approaches is in our assumption of complete ionization, while Storch et al. (2014) used a partially ionized atmosphere. Surprisingly, the effect on the pulsed fraction is quite substantial, despite the fact that the fraction of neutral H atoms for the typical values of the temperature and magnetic field relevant to PSR B0943+10 is well below 0.1% (Potekhin & Chabrier 2004).

Given the uncertainties in parameters (${T}_{\mathrm{eff}}$, B) and in the physical conditions of the atmosphere, one may tentatively conclude that magnetic beaming in an atmosphere on top of the heated polar cap can reproduce the large observed pulsed fraction, and its dependence on energy, for values of the two geometrical angles ξ and χ in the range derived by radio observations.

Finally, having achieved a reasonable characterization of the X-ray components in the B- and Q-modes, we turned our attention to see whether our extensive observations could tell us anything about the evolution of these characteristics within the modes. A noticeable change in the shape of the pulsar radio profile and in the drifting feature occurs during the course of the B-mode (Backus et al. 2011; Bilous et al. 2014). The work by Backus et al. (2011) revealed that, based on the carousel model, the change in drift could be produced by a 10% increase in the average number of sub-beams and a 16% increase in the circulation time of the carousel. They speculated that under the PSG model the increase in circulation time should be related to a change in temperature of about 1.4% across the B-mode. Such a small change is clearly consistent with our X-ray data, but unfortunately cannot be detected with the current instrumentation.

In Section 3.4 we presented marginal evidence for an increase in the pulsed fraction during the B-mode. We note that, if this occurs on the same timescale as the evolution of the radio properties, our analysis underestimates its significance since the optimal boundary, with more plentiful X-ray counts, would have been, say, 1.5 hr rather than the selected 3 hr. The failure to detect pulsations in the B-mode during the 2011 observations, in which the later stages of the B-mode were poorly covered (Hermsen et al. 2013), might be taken as a further indication in favor of an evolution of the pulsed flux. If this effect is confirmed with more X-ray data, its study can provide interesting clues for the understanding of the processes responsible for the mode-switching behavior of PSR B0943+10.