COORDINATED X-RAY, ULTRAVIOLET, OPTICAL, AND RADIO OBSERVATIONS OF THE PSR J1023+0038 SYSTEM IN A LOW-MASS X-RAY BINARY STATE

PSR J1023+0038 (also known as AY Sextantis or FIRST J102347.6+003841) is a 1.7 ms pulsar in a 4.75 hr binary orbit around a bloated ∼0.2 ${M}_{\odot }$ main-sequence-like companion star. This system is notable in that it was the first to exhibit compelling evidence for the transition process between an accretion disk-dominated low-mass X-ray binary (LMXB)-like state and a disk-free radio pulsar state. Optical observations revealed an accretion disk in the system in 2001 (Bond et al. 2002; Szkody et al. 2003; Wang et al. 2009) that seemed to be absent in early 2003 (Thorstensen & Armstrong 2005) and at the time of the radio pulsar discovery (Archibald et al. 2009). During its disk episode in 2001, the binary appeared optically blue and bright in addition to showing strong, double-peaked H and He emission lines, which are commonly seen in LMXBs, and are a typical feature of accretion disks. This suggests that the companion star overflowed its Roche lobe, forming a disk surrounding the millisecond pulsar (MSP). In contrast, spectra obtained starting in early 2003 showed a substantially lower optical flux and a typical G-type spectrum, implying the absence of a substantial accretion disk.

A recent torrent of developments have revealed two close analogs to PSR J1023+0038. In 2013 March, PSR J1824– 2452I in the globular cluster M28 was seen to switch between rotation-powered (radio) and high-luminosity accretion-powered (X-ray) pulsations (Papitto et al. 2013), thereby strenghtening the long-suspected evolutionary link between LMXBs and "recycled" pulsars (Alpar et al. 1982). In addition, re-examination of archival optical and X-ray data of the X-ray binary XSS J12270–4859 revealed that only a few months prior (in 2012 November/December) its accretion disk had disappeared and the 1.7 ms radio pulsar became active (Bassa et al. 2014; Bogdanov et al. 2014b; Roy et al. 2014). See Linares (2014) for an overview of the X-ray states of these and analogous objects.

Unexpectedly, on 2013 June 23, no radio pulsations were detected from PSR J1023+0038, and none have been detected up to the publication of this article despite an intensified monitoring campaign and higher-radio-frequency observations with greater sensitivity. Subsequent X-ray, optical, and γ-ray investigation revealed that PSR J1023+0038 has undergone another transformation to an accretion disk-dominated state (Halpern et al. 2013; Patruno et al. 2014; Stappers et al. 2014; Tendulkar et al. 2014). This transition was accompanied by an extraordinary five-fold increase in high-energy γ-ray luminosity as seen by the Fermi Large Area Telescope (Stappers et al. 2014). The X-ray emission observed by Swift XRT revealed a significant change in behavior compared to the disk-free state as well: in addition to an increase in flux by over an order of magnitude to ∼10^33–34 erg s⁻¹, the orbital-phase-dependent modulations had disappeared, replaced by aperiodic variability with rapid drops to a low flux level (Patruno et al. 2014). Follow-up hard X-ray observations with NuSTAR revealed strong flares reaching up to $\approx 1.2\times {10}^{34}$ erg s⁻¹ (3–79 keV), as well as the same peculiar drops in flux (Tendulkar et al. 2014). This behavior bears close resemblance to that observed in PSR J1824–2452I (Linares et al. 2014) and XSS J12270–4859 (de Martino et al. 2010, 2013) in their low-luminosity LMXB states.

In Archibald et al. (2015), we established that despite PSR J1023+0038 existing at a luminosity level typical of quiescent LMXBs, channeled accretion onto the neutron star surface appears to be occuring on a regular basis as demonstrated by the detection of coherent X-ray pulsations. Herein, we build upon this finding with a more in-depth analysis of the same XMM-Newton data set, which is augmented by an array of coordinated multi-wavelength observations of PSR J1023+0038. This study provides the most complete observational picture to date of the LMXB state of the PSR J1023+0038 system and, arguably, transition MSP systems in general. The work is organized as follows. In Section 2, we summarize the observations and data analysis. In Section 3, we focus on the X-ray variability, while in Section 4 we focus on its statistical properties. We discuss the X-ray pulsations in Section 5 and summarize the spectroscopic analysis of the X-ray emission in Section 6. We describe the results of the optical photomety and spectroscopy in Sections 7 and 8. We present the results of our radio imaging and timing observations in Sections 9 and 10. Finally, in Section 11 we provide a discussion and offer conclusions in Section 12.

Tables 1 and 2 summarize all X-ray, UV/optical, and radio observations presented in this paper, while Figure 1 shows a visual timeline of this multi-wavelength observing campaign, arranged around two XMM-Newton observations conducted in 2013 November and 2014 June.

Zoom In Zoom Out Reset image size

2.1. X-Ray Observations

2.2. Optical/UV Photometry

2.3. Optical Spectroscopy

2.4. Radio Timing Observations

2.5. Radio Imaging Observations

Figures 2 and 3 show the total exposure-corrected, background-subtracted XMM-Newton EPIC X-ray light curves in the 0.3–10 keV band, obtained by combining the data from the pn, MOS1, and MOS2 during the periods when all three telescopes acquired data simultaneously. The large-amplitude variability is obvious and is reminiscent of the behavior seen in the Swift XRT (Patruno et al. 2014) and NuSTAR (Tendulkar et al. 2014) observations of PSR J1023+0038 from 2013 October. During both XMM-Newton obervations, most of the time is spent in the "high" mode,¹⁰ with a typical total EPIC count rate of ≈7 counts s⁻¹ (0.3–10 keV). The emission drops out unpredictably to the "low" mode with count rate ≈1 counts s⁻¹ with very rapid ingress and egress (∼10–30 s).

Zoom In Zoom Out Reset image size

Sporadic, intense X-ray flares (labeled by the Roman numerals in Figures 2 and 3), reaching up to ≈60 counts s⁻¹, occur on average every few hours, and exhibit diverse morphologies and durations. Some flares last less than a minute (e.g., II), while others (e.g., III, IX, and XI) last up to ≈45 minutes. The long-duration flares exhibit a great deal of intricate structure. Specifically, throughout flares III and XI, there are a number of rapid drops in count rate that occasionally reach down to the levels of the low flux mode. In both XMM-Newton observations, the reccurence time between X-ray flares is comparable (∼20 ks), with six prominent flares (which we define as those with peak rates exceding 25 counts s⁻¹) in each observation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is interesting to note that the brighter flares are usually (but not always) preceded by a low mode instance (see, e.g., flares I, IV, V, VI, VII, VIII, X, and XI). This behavior could be an indication of a physical process that depends on the accumulation of a reservoir of energy (during the low mode) that is rapidly released (in the flare). However, due to the limited number of strong flares in the present data combined with the fairly frequent occurrence of low mode intervals, it is not possible to determine with certainty if this is merely a chance coincidence or if the flares and their corresponding pre-flare low mode are truly associated. The absence of low modes immediately preceding some flares argues against a causal connection because it implies that the occurence of a low mode interval is not required to initiate a flare.

In the disk-free, radio pulsar state, PSR J1023+0038 exhibited pronounced orbital-phase-dependent X-ray flux modulations (Archibald et al. 2010; Bogdanov et al. 2011). Although in the light curves shown in Figures 2 and 3 there is no obvious periodic behavior (see Section 4), in principle, orbital modulation may still be present in one of the X-ray flux modes. To investigate this possibility, we separately extracted the emission from the high, low, and flare modes and folded them at the orbital period. We do not find evidence for orbital modulation of the X-ray brightness during any of the three modes.

Li et al. (2014a) have reported evidence for periodicity at 3130 s in the NuSTAR data set of PSR J1023+0038. We have examined both XMM-Newton observations in search of a similar signal but find no evidence for it. In addition, we are not able to reproduce this result using the NuSTAR data presented in Tendulkar et al. (2014).

We performed a statistical analysis of the low, high, and flare X-ray flux modes to understand their temporal properties and identify any correlations between them. We followed the same methodology as the analysis performed on the NuSTAR data of PSR J1023+0038 in Tendulkar et al. (2014). We use the combined background-subtracted and exposure corrected EPIC 0.3–10 keV data binned with 10 s bins shown in Figures 2 and 3. The 10 s bins were chosen to ensure sufficient photon statistics in the low mode bins and in the mode transitions.

Figure 5 shows the distribution of count rates throughout the 2013 November (solid black line) and 2014 June (dashed red line) observations. The bi-modality in the distribution is immediately apparent, with the low and high modes clearly seen as peaks at ≈1 counts s⁻¹ and ≈7 counts s⁻¹, respectively, for both epochs. In both instances, the distribution of count rates for the low and high mode are consistent with what is expected from a Poisson distribution. The minimum of the distribution between the two states occurs at approximately 3.1 counts s⁻¹.

Zoom In Zoom Out Reset image size

To cleanly differentiate between the two modes and the transitions between them, we define a "gray area" (hatched region in Figure 5) ranging from 2.1 to 4.1 counts s⁻¹. These limits were chosen to be symmetrical around the minimum of the distribution, though because of the inherent count-rate errors in each light curve bin, the exact values of these thresholds is not critical to the analysis presented below.¹¹ In particular, varying the thresholds by ±0.3 s⁻¹ results in a difference of less than 1% in the number of transitions. Figure 5 shows the count rate levels used to define the low and high states. In the analysis that follows, the flare mode intervals were removed from the light curves by excising intervals with rates in excess of 11 counts s⁻¹.

We use a bi-stable comparator¹² to define the transitions between the modes as follows. A low→high transition is defined when the count rate crosses from the "low" zone to the "high" zone shown in Figure 5 and reverse for the high→low transition. If the light curve varies from the low region to the gray area and returns back to the low region, the light curve mode is maintained as low and a transition is not registered. Similarly, a variation from the high zone to the gray area and back to the high zone is maintained as a high mode. Using this comparator algorithm, we counted 237 low→high and high→low transitions in the 2013 November data and 172 low-high transitions and 171 high-low transitions in the 2014 June data. We measured the durations of the low and high modes as the number of light curve bins between the end of the previous transition and the beginning of the next transition. The duration of the transition itself was estimated as the number of 10 s light curve bins spent in the gray area.

Figure 6 shows the distribution of durations of the low (dashed red line) and high (blue line) mode based on the combined 2013 November and 2014 June data sets. The histogram is truncated to show only durations shorter than 1000 s although the longest continuous high mode interval was 6910 s in 2014 June and the longest low mode lasted 1930 s, also in 2014 June. We find that the low mode exists for $\approx 22$% of the total observation time during 2013 November 10–12 and $\approx 21$% during 2014 June 10–11 observations. As shown in Figure 6, on average the high mode tends to persist for significantly longer durations than the low mode. The distribution of low mode durations can be very approximately modeled as a power law with index −1.1 and −1.3, for the 2013 November and 2014 June observations, respectively. The separations of the centroids of consecutive low modes follow a log-normal distribution with mean separations of ∼441 and ∼403 s and standard deviations of 0.93 and 0.76 for the two data sets, respectively.

Zoom In Zoom Out Reset image size

We also looked for correlations or patterns between durations of consecutive modes. There is no evidence for correlated behavior between the two modes, indicating a stochastic nature of the underlying physical mechanism responsible for the switches between the two. However, we do note that the modes tend to last longer during 2014 June as compared to the durations in 2013 November, which suggests some long-term variation in behavior.

We investigated the distribution of the times of the low-high and high-low transitions for 2013 November and 2014 June as well. We find that 57% of the high→low transitions in 2013 November and 54% of those in 2014 June lasted <10 s. In contrast, only 40% of the low→high transitions in 2013 November and 39% of those in 2014 June were <10 s. Thus, based on this estimate, the high→low transitions may be faster than the low→high transitions.

As reported in Archibald et al. (2015), the XMM-Newton data revealed coherent modulations at the pulsar spin period in both the 2013 November and 2014 June observations. Intriguingly, the pulsations occur just in the high flux mode, with a pulsed fraction of 8% (0.3–10 keV). The profile is clearly double-peaked, with a ∼180° separation in rotational phase, with each peak presumably corresponding to emission from one of the diametrically opposite magnetic polar caps of the neutron star.

Here we examine the energy dependence of the high mode X-ray pulse profiles as well as any spectral variations as a function of spin phase. The method for extracting the pulse profiles is given in Archibald et al. (2015). The top right panel of Figure 7 shows the high mode pulse profile in four energy bands (0.3–0.7, 0.7–1.5, 1.5–3, and 3–10 keV) obtained by aligning the total profiles from the two individual observations via cross-correlation (see Archibald et al. 2015, for details). The hardness ratio of the high mode exhibits a subtle spectral softening in the trailing edge of the stronger pulse. For reference, we also show the flare and low mode data folded using the same ephemeris, where no statistically significant pulsations are seen. The strongest, albeit marginal (with a single-trial significance of only ∼2σ) signal is found above 3 keV for the flare mode, although it shows only single-peaked modulation.

Zoom In Zoom Out Reset image size

As shown in Figure 8, the spectrum of the pulsations is similar to the spectrum of the source generally. This suggests that the pulsed and most of the unpulsed radiation are produced by the same process. It is possible that the pulsations decline above 5 keV but this cannot be conclusively established due to the limited photon statistics above this energy.

Zoom In Zoom Out Reset image size

Previous studies of PSR J1023+0038 have established that the X-ray emission in the LMXB state is well-described by an absorbed power-law with spectral photon index of $\rm{\Gamma }\approx 1.7$ with no substantial spectral changes despite the over one order of magnitude variation in flux (Coti Zelati et al. 2014; Patruno et al. 2014; Takata et al. 2014; Tendulkar et al. 2014). The large photon harvest of the XMM-Newton data presented here allows us to identify any subtle differences in spectral properties between the three distinct flux modes: low, high, and flare. The corresponding spectra were obtained using selections based on the count rate ranges for the low and high modes determined in Section 4, while for the flares we choose a count rate cut of ≥15 counts s⁻¹ to minimize contamination from the high mode. Due to the relatively large uncertainties in the spectral calibration of pn timing mode data (see, e.g., Archibald et al. 2010; Bogdanov 2013), we only consider the MOS1/2 spectra in this analysis.

Table 4 summarizes the results of the power-law spectral fits for the three flux modes fitted in XSPEC (Arnaud 1996). We conducted the analysis in two ways: one with each mode fitted separately and the other with all three modes fitted jointly with a tied value of the hydrogen column density along the line of sight (N_H). In both instances, the fits suggest that the value of N_H is consistent among the three with no appreciable enhancement in the intervening absorbing column within the binary during any of the modes.

Table 4.  Results of Absorbed Power-law Fits of the XMM-Newton Spectra of PSR J1023+0038

	NH	Γ	LXb	${\chi }^{2}$/dof
Mode	(1020 cm−2)	⋯	(1033 erg s−1)
2013 Nov 10–12 (ObsID 0720030101)
Flare	2.7(9)	1.65(4)	10.8(2)	0.98/321
High	3.1(2)	1.71(1)	3.17(2)	1.23/902
Low	2.6(11)	1.80(5)	0.54(1)	1.02/281
Flarec	3.1(2)	1.66(2)	10.9(2)	1.14/1506
High	⋯	1.71(1)	3.17(2)	⋯
Low	⋯	1.82(3)	0.54(1)	⋯
2014 Jun 10–11 (ObsID 0742610101)
Flare	3.6(8)	1.76(3)	9.6(2)	1.02/373
High	3.1(2)	1.75(1)	3.06(2)	1.31/870
Low	2.4(12)	1.77(6)	0.45(1)	0.96/206
Flarec	3.1(2)	1.74(2)	9.6(1)	1.18/1452
High	⋯	1.75(1)	3.06(2)	⋯
Low	⋯	1.80(4)	0.46(1)	⋯

Notes.

^aThe numbers in parentheses show the 90% confidence level uncertainties in the last digit of the quoted best-fit values.

^bLuminosity in the 0.3–10 keV band assuming the parallax distance of 1.37 kpc (Deller et al. 2012). ^cFits performed jointly on the spectra of all three modes with a tied value of N_H.

Download table as:  ASCIITypeset image

For the flare and low modes, the emission is well-described by an absorbed power-law in both observations. In contrast, the fits to the high mode spectra yield null hypothesis probabilities of 3 × 10⁻⁶ and 2 × 10⁻⁹ for the 2013 November and 2014 June data, respectively. This means that a simple absorbed power-law model alone does not provide an adequate representation of the high mode spectrum, which is also evident from the broad-band residuals shown in Figure 9. One likely explanation is the presence of a faint thermal (e.g., neutron star atmosphere) component since the X-ray pulsations imply accretion onto the neutron star, which may result in at least some superficial heating of the surface. To account for this emission we consider two models: nsa, a passive, weakly magnetized neutron star hydrogen atmosphere (Zavlin et al. 1996) and, zamp, a hydrogen atmosphere accreting at a low rate (Zampieri et al. 1995). For the latter, the spectrum is defined in terms of the observed accretion luminosity expressed in terms of the Eddington rate, assuming radiation from the entire surface of a neutron star with ${M}_{\mathrm{NS}}=1.4$ ${M}_{\odot }$ and R_NS = 12.4 km. The best fit parameters of the thermal component are consistent between the 2013 November and 2014 June observations (Table 5). The implied effective emission radii of the nsa model are smaller than the whole neutron star, consistent with emission from hot spots. The implied contribution of thermal radiation to the total luminosity (3%–9%) is consistent with the pulsed fraction of the X-ray pulsations in the high mode. The implied luminosity from the zamp model is $2.5\times {10}^{-5}{L}_{\mathrm{edd}}$, which translates to a thermal luminosity of ∼6 × 10³³ erg s⁻¹, much greater than the implied thermal fraction. This arises due to the fact that the model considers emission from the entire surface, whereas the pulsations from PSR J1023+0038 indicate emission from a small portion of the surface. Although not strictly correct, adjusting the size of the emitting area to that obtained with the nsa model produces a value consistent with the thermal luminosity deduced from the spectral fits. While both composite models result in a reduction of the broad residuals relative to the pure power-law model, for both sets of observations the null hypothesis probabilities (7 × 10⁻³ and 1 × 10⁻³, respectively) indicate fits that are formally not acceptable. It may be that the flare and low state have similar residuals but due to the lack of photon statistics they are not as apparent.

Zoom In Zoom Out Reset image size

Table 5.  Multi-component Model Fits of the High Mode XMM-Newton Spectra of PSR J1023+0038

	NH	Γ	Teff	Reff	${L}^{\infty }/{L}_{\mathrm{Edd}}$	LX,PLb	LXb	${\chi }^{2}$/dof
Modela	(1020 cm−2)		(106 K)	(km)		(1033 erg s−1)	(1033 erg s−1)
2013 Nov 10–12 (ObsID 0720030101)
powerlaw + nsac	2.5(6)	1.59(3)	1.49(21)	2.9 + 2.6−1.9	⋯	3.00(4)	3.18(3)	1.12/900
powerlaw + zamp	2.5(5)	1.60(3)	⋯	⋯	−4.6(2)	3.02(4)	3.18(3)	1.12/900
2014 Jun 10–11 (ObsID 0742610101)
powerlaw + nsa	2.3(6)	1.60(3)	1.45(18)	${3.3}_{-1.9}^{+2.7}$	⋯	2.85(4)	3.06(3)	1.15/868
powerlaw + zamp	2.2(5)	1.61(3)	⋯	⋯	−4.6(2)	2.87(4)	3.06(3)	1.15/868

Notes.

^aAll uncertainties quoted are at a 90% confidence level. ^bLuminosity in the 0.3–10 keV band assuming the parallax distance of 1.37 kpc (Deller et al. 2012). ^cFor the nsa model a neutron star with M = 1.4 ${M}_{\odot }$ and R = 12 km is assumed.

Download table as:  ASCIITypeset image

A possible explanation for the substantial residuals in the high mode fits could be the presence of a complicated spectrum produced by three or more distinct emission components. However, adding components commonly used in modeling of LMXB spectra such as thermal Comptonization or bremsstrahlung does not produce acceptable fits either.

An alternative interpretation involves small fluctuations of the power-law photon index during the high mode, which would result in a mean spectrum that is not well represented by a power-law with a single photon index. To investigate this possibility we have examined the ratio of counts in the 1.5–10 keV and 0.3–1.5 keV bands versus count rate using the combined pn and MOS1/2 data. We consider three count rate ranges within the high mode with bounds 4.7, 6.2, 7.8, and 9.3 counts s⁻¹, chosen to minimize the contamination from the low and flare modes (see Figure 5). For the 2013 November data, this hardness ratios are 0.666 ± 0.003, 0.678 ± 0.003, and 0.804 ± 0.005 at the lower end, peak, and upper end of the count rate distribution of the high mode, respectively. For the 2014 June data, the values are 0.637 ± 0.002, 0.606 ± 0.002 and 0.685 ± 0.005, respectively. In both instances, the difference in hardness as a function of count rate is highly significant, which is an indication of spectral variability within the high mode. Moreover, for a given count rate the hardness ratio is not the same between the two observations with the spectrum being generally softer in the 2014 June data (as evident from the best-fit power-law indices in Table 4 and confirmed by the values of the hardness ratio). Therefore, fluctuations in the power-law index can explain the poor quality of the spectral fits.

The flux in the flare, high, and low modes appears to decrease by 11%, 3.5%, and 16% between 2013 November and 2014 June. We have investigated whether this could be due to the choice of count rate cuts for the three flux modes. Using more conservative cuts produces similar results for the flux and photon index. Thus, the long-term flux changes of the modes appear to be intrinsic to the system.

We note that for the case of the flare spectra the quoted luminositites represent time-averaged values. The peak luminosity, observed in flare IX, reaches a luminosity of ≈3 × 10³⁴ erg s⁻¹ (0.3–10 keV), assuming the same spectral shape as the total flare spectrum.

As noted by Halpern et al. (2013), in the optical the PSR J1023+0038 binary is ∼1 magnitude brighter during the LMXB state than in the radio pulsar state. From the XMM-Newton OM and MDM photometric light curves (Figures 10 and 11) orbital modulation similar to the heating light curve observed during the radio pulsar state (Woudt et al. 2004; Thorstensen & Armstrong 2005) are also evident. However, the variations now appear more sinusoidal and symmetric about the peak. The optical brightness varies between $B\approx 17.5$ and $B\approx 16.7$, although orbit-to-orbit evolution of the overall brightness is apparent. For comparison, in the radio pulsar state the B filter magnitude varies between 18.4 and 17.8 (Thorstensen & Armstrong 2005; Homer et al. 2006).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A wide variety of flaring behavior is also seen, which is associated with the accretion disk state. Some of these flares are very rapid, lasting less than 1 minute, with a rise time that is unresolved at our 10 s (OM) and 13 s (MDM) cadences. A good example of this phenomenon can be seen in the MDM light curve on December 26.51 UT. At the other extreme is a continuous episode of strong flaring on December 28 that lasts 3.6 hr, or 3/4 of the binary orbit. Low-level flaring is present almost continuously for the entire 2013 November XMM-Newton observation and the MDM 2013 January 3 observation. We do not see any optical moding behavior that resembles the X-ray variability. In all cases, the optical variability can be characterized as positive flares superposed on the otherwise smooth heating light curve.

The B filter photometric data obtained with the XMM-Newton OM further reveal that the brightest optical flares closely match the prominent X-ray flares. Therefore, the flares in PSR J1023+0038 appear to be broad-band phenomena that span from the optical up to at least hard X-rays (as seen with NuSTAR; Tendulkar et al. 2014). Of particular note are flares III and XI, which reach peak B filter magnitudes of 15.9 and 15.5, respectively. For some, especially flares I and V, the flux peak occurs at significantly later times than the X-ray one and the optical tail of the flare is much longer (see Figure 4). For the brief but intense X-ray flare labeled as II, there does not appear to be a corresponding optical event, while multiple fainter rapid X-ray flares (with X-ray count rates below 20 s⁻¹) seen in Figures 2 and 3 have clear optical counterparts.

To identify any time lags between the X-ray and optical flares, we cross-correlated the XMM-Newton X-ray and optical time series binned at 10 s resolution. We first fitted a sine function to the OM observation and subtracted it from the data to remove the orbital modulation. The best-fit sine function results in a period of 4.754 hr, consistent with the orbital period of the system. As expected, the cross-correlation between the X-ray and detrended optical data sets reveals a strong correlation, which is dominated by the flares. The time lags among the flares range between −260 and +160 s, where a positive lag indicates that the X-ray flare precedes the optical one. The strongest correlation is found for the short flares I (with time lag +130 s) in 2013 November and X (with lag 0 s) in 2014 June. For flares III, VIII, IX, and XI the optical leads the X-ray (by −160, −260, −100, −30 s, respectively) although the correlation is weak. The time lags of individual flares are unlikely to be correlated with the orbit due to the variation of the emission site across the orbit with respect to the observer because the compactness of the binary (∼4–5 light s). Indeed, we find no clear trend of the time-lag changes with the orbital phase. Finally, we investigated any correlation of the low level flux variation between the X-ray and optical by removing all significant flares. The resulting correlations are very weak, consistent with the apparent lack of flux moding in the optical photometric data.

Another way to reveal low-level optical flux variations correlated with the X-ray mode switching is to consider the optical emission averaged over many low and high X-ray mode intervals. To this end, we first removed the sinusoidal variations of the optical time series. For each X-ray low mode we then computed averaged B filter light curves for specific durations (100, 200, and 300 s) and compared them against the high mode immediately before and after using the same duration. The 2013 November data show that the source is on average ∼0.1 mag brighter during the X-ray low modes relative to the high modes that occur immediately before or after. However, the 2014 June data does not show the same behavior.

The UV emission observed with the Swift UVOT (see Figure 12) exhibits rapid variability as well. A moderately bright flare that occurs around MJD 56606.83 has an obvious UV counterpart in the Swift UVOT data (as seen in the left-most panels of Figure 12). The UVOT UVW1 filter data from 2013 November and the UVW2 filter data from 2014 June cover time periods when multiple low flux intervals are seen in the XMM-Newton data. There are no obvious UVW1 brightness variations commensurate with the nearly order of magnitude X-ray flux changes associated with the mode transitions. The 2014 June UVW2 data exhibits a brightness increase that coincides with a low-high X-ray mode transition, although due to the brief UVOT exposure it is not clear if this flux variation occurs consistently for mode transitions.

Zoom In Zoom Out Reset image size

In Figure 13, the average of the two individual VLT spectra is plotted. The spectrum is dominated by a blue continuum and strong, double-peaked emission lines of H and He. The spectrum appears similar to other optical spectra of the PSR J1023+0038 binary presented in the literature, both during the active period in 2000/2001 (Bond et al. 2002; Wang et al. 2009) and 2013/2014 (Halpern et al. 2013; Coti Zelati et al. 2014; Takata et al. 2014). The optical spectra confirm that the accretion disk was present shortly after the end of the 2013 November X-ray observation. As the optical spectra presented in the literature since the 2013 June transition all show the presence of an accretion disk, it is highly probable that it has been present throughout the 2013 November and 2014 June X-ray observations of PSR J1023+0038. A detailed analysis of these and a series of other spectra acquired since the state transformation of PSR J1023+0038 will be presented in a subsequent paper.

Zoom In Zoom Out Reset image size

As reported in Stappers et al. (2014), the 2013 June reappearance of the accretion disk in PSR J1023+0038 was accompanied by cessation of its bright radio pulsations. To verify this using the 2013 November and 2014 June observations, the radio data were all folded with the ephemeris derived prior to the pulsar disappearance in 2013 June. In both cases sub-integration times of 10 s were used with 1600 frequency channels of with 0.25 MHz for the LT and 512 frequency channels of width 0.3125 MHz for the WSRT data. No radio pulsations were detected around the pulsar period.

Conceivably, if the rotation-powered pulsar emission mechanism is actually still active, the radio pulsations may be highly intermittent such that they may only be present in one of the three distinct X-ray flux modes. In such a scenario, the radio emission is either quenched or severely obscured in the other modes and is only able to "break through" during one mode. To investigate this possibility, we divided the radio observations into time segments that match the times of the high, low, and flare modes observed during the simultaneous XMM-Newton observations. For this purpose, we defined an X-ray low mode when the count rate is less than 3 counts s⁻¹ and the duration of the low mode is more than 40 s. This resulted in 70 and 35 low mode intervals overlapping with the LT and the WSRT data, respectively. The data in these sections were then corrected for the new, local ephemerides derived from the X-ray data (from Archibald et al. 2015) and subsequently a search in period and dispersion measure (DM) was undertaken for each section. The period search range was 0.01 μs and the DM search range was 3 cm⁻³ pc. These values were chosen to allow for significant changes in the orbital parameters and/or in the electron column along the line of sight due to material being lost from the companion, or due to the accretion disk. To a signal-to-noise limit of 6 we find no evidence of pulsed emission during any of the low mode intervals in either the LT or WSRT data. We note that these data span a wide range of orbital phases when the pulsar has previously been detected with S/N of many tens in similar integration times. We performed the same analysis on flares that occur during the radio observation with no significant detection of pulsations either. Finally we considered a long sequence of high mode X-ray emission (MJD 56607.348–56607.385) overlapping with the LT data and also find no significant pulsed signal.

The two-hour long VLA observation made at 5 and 7 GHz at the beginning of the 2013 November campaign shows a strong detection, with a flat spectrum ($\alpha =0.09\pm 0.18$, where flux density is proportional to ${\nu }^{\alpha }$) and flux density which varies rapidly between 60 and 280 μJy (see left panel of Figure 14). In the 5 GHz e-EVN VLBI observations made shortly after the 2013 November observations we obtain a marginal detection (3.8σ peak with flux density 50 μJy at the precise position of the pulsar predicted by the astrometric solution of Deller et al. 2012); the much sparser uv coverage and the faintness of the source precluded any detection of variability on a sub-observation timescale. In the shorter 2–4 GHz VLA observation made 20 hr later, the source is not detected in the combined dataset, with a 3σ upper limit of 30 μJy. When a light curve is made with a resolution of 4 minutes, two tentative, 3σ detections with peak flux density ∼60 μJy are made (right panel of Figure 14). Finally, PSR J1023+0038 is not detected in the 140 MHz LOFAR observations, although the upper limit of several mJy is not constraining (we note however that, when visible, the pulsar has a period-averaged flux density of ∼40 mJy at 150 MHz; V. Kondratiev et al. 2015, in preparation). Deller et al. (2014) show that the variable, flat-spectrum emission seen here persists over many observations made in a six-month period after the 2013 November observations. The observed radio emission is inconsistent with a radio pulsar origin, which should be very steep ($\alpha \sim -2.8$; Archibald et al. 2009) and hence much brighter than the observed emission below 4 GHz, while being visible up to ∼6 GHz.

Zoom In Zoom Out Reset image size

Only the last 16 minutes of the VLA data from 2013 November 10 overlap with the XMM-Newton EPIC MOS1 exposure.¹³ During most of this period the X-ray emission is in the low flux mode, while the radio emission exhibits a rise and fall in flux; the overlap is too brief to establish whether any correlation (potentially with a time lag due to different production sites for the radio and X-ray emission) exists. Although the XMM-Newton OM data cover most of the VLA observations, no discernable relation between the radio and optical variability is present. For the 2013 November 11 radio observation, the faintness of the radio emission makes it impossible to identify any relation between the radio and X-ray emission.

Although it is thus difficult to establish whether there is any correlation between the radio and X-ray light curves, we can still examine whether there are any noticeable differences in the radio emission during the different X-ray modes. As discussed in Section 11.1, this has significant implications for whether the radio pulsar could still be active at any time during the current LMXB state, since the activation of the radio pulsar would lead to the presence of additional steep-spectrum radio emission. A significant portion of the 16 minutes overlap between the MOS1 and VLA observations on November 10 is in the low mode, and when we image the VLA data from 18:00:00 UT to 18:06:30 UT (during this mode) separately, we obtain a flux density of 180 ± 15 μJy with a spectral index of 0.4 ± 0.6: consistent with the average value for the whole observation. During the 2–4 GHz observation on 2013 November 11, multiple examples of the low mode are visible, generally for a period of a few minutes at a time. In a radio light curve with 1 minute resolution, no peaks >2.5σ (∼100 μJy) were visible. In either of these cases (the 6.5 minutes low mode image in the November 10 observation and the 1 minute time slices in the November 11 observation) the radio pulsar emission, if unchanged compared to that reported in Archibald et al. (2009), would have been clearly visible. We discuss the implications in Section 11.1.

Our unprecedented multi-wavelength campaign of PSR J1023+0038, conducted in 2013 November and 2014 June, shows fascinating phenomenology that provides fresh insight into the properties of PSR J1023+0038 and analogous systems. Based on the wealth of data we have accumulated, we are able to establish the following.

In its LMXB state, PSR J1023+0038 resides in a luminosity range of ${L}_{X}\sim {10}^{32-34}$ erg s⁻¹ (0.3–10 keV). The majority of the time is spent at a $3\times {10}^{33}$ erg s⁻¹ level, during which coherent X-ray pulsations are observed. This implies that active accretion onto the stellar surface takes place (Archibald et al. 2015). During the low flux mode intervals ($\sim 5\times {10}^{32}$ erg s⁻¹), no pulsed emission is seen, possibly because the accretion flow is unable to reach the star. As shown in Archibald et al. (2015), there is no evidence for pulsations in the low flux mode, with a 95% confidence upper limit on the pulsed fraction of ∼2.4%. This excludes the possibility that the low mode simply has the same amplitude of pulsations but with scaled-down luminosity. Instead, it appears that the coherent X-ray pulsations are either completely absent or have a comparatively very low amplitude.

The low mode 0.3–10 keV luminosity is ∼5 times greater than the average luminosity observed in the radio pulsar state in the same energy band (9 × 10³¹ erg s⁻¹). In the radio pulsar state, PSR J1023+0038 showed evidence for broad double-peaked pulsations at the spin period at a 4.5σ with a ∼9% pulsed fraction for 0.3–10 kev (Archibald et al. 2010), in addition to the dominant shock emission component. Such pulsations are observed in multiple rotation-powered MSPs and are believed to be produced by heating of the magnetic polar caps by a return flow of relativistic particles from the pulsar magnetosphere. If the X-ray pulsations seen in the radio state are present in the low mode they would appear with a pulsed fraction of ∼2%, below the upper limit we have derived. Therefore, we cannot rule that the same X-ray pulsations seen in the radio state are present in the low mode of the LMXB state.

Archibald et al. (2015) demonstrated that the pulsations appear to be completely absent during the flares as well. The pulsed fraction upper limit of 1.5% corresponds to a pulsed flux of 0.21 counts s⁻¹, lower than the pulsed flux in the high flux mode (0.33 counts s⁻¹). This indicates that during the flares the pulsations are strongly supressed. Therefore, the flares are not simply added flux on top of the high mode emission as the pulsations would still be easily detectable. Instead, the high mode appears to be absent during the flares.

The multiple intense X-ray flares exhibit varied morphologies and temporal behavior. Flaring/burst activity in accreting systems might be from coronal flares arising due to magnetic reconnection events in the accretion disk (e.g., Galeev et al. 1979). Alternatively, an instability in the inner disk may cause a large influx of acreting material toward the compact object (see Section 11.4). Of particular interest are the long flares interspersed with rapid drops in flux (e.g., III an IX in Figure 4) down to X-ray flux levels comparable to the low mode. This suggests that at least for these flares the emission is not superposed on the high mode flux, which offers an additional line of evidence for the lack of high mode flux in the flares.

There is no indication for the orbital-phase-dependent X-ray modulations that were seen in the radio pulsar state even if all intervals of a single flux mode are folded at the binary period. The non-thermal spectrum in each of the three flux modes is much softer compared to that seen in the radio pulsar state, which had Γ = 1 − 1.2 (Archibald et al. 2010; Bogdanov et al. 2011). This implies that the intra-binary shock near the face of the companion, responsible for the X-ray emission in the disk-free state, is either absent or completely overwhelmed in the accreting state. The change in the optical heating light curve from an asymmetric to a symmetric, more sinusoidal shape favors the absence of this shock emission in the LMXB state. This difference can be explained if the heating by the pulsar wind was asymmetric due to channeling of the pulsar wind by the companion's magnetic field (see, e.g., Li et al. 2014b), while the heating in the LMXB state is by X-rays from the inner disk, which would produce a more symmetric light curve as the source of heating is more point-like. In addition, even if the rotation-powered pulsar mechanism is still active in the LMXB state, the accretion flow may intercept a substantial fraction of the pulsar wind directed toward the companion (see, e.g., Takata et al. 2014).

The sudden drops to a low flux mode appear to be of an entirely different nature than the dips in the X-ray dipper variety of LMXBs (e.g., White & Swank 1982; Smale et al. 1988; Balucinska-Church et al. 2011) as they do not show spectral changes in either soft X-rays (as shown in Section 6) or hard X-rays (see Tendulkar et al. 2014). In general, based on the results of the X-ray spectral fits presented in Section 6, any interpretation that evokes obscuration/absorption is unlikely. Since there are negligible spectral changes, especially the value of N_H (see Tables 3 and 4), between the high and low modes, there is no indication for an enhancement in the amount of intervening material during the low mode. There is also no evidence of an X-ray luminosity dependence on the duration and frequency of flares and low flux mode intervals or any correlation between the separation between (and duration of) dips or flares.

Although the extensive data set presented here reveals a wide variety of interesting flaring phenomenology there is evidence that it does not cover the whole range of possible behavior of PSR J1023+0038 in the LMXB state. In particular, in the 100 ks NuSTAR observation of PSR J1023+0038 from 2013 October, the system undergoes a prolonged episode of flaring activity that lasts for ∼10 hr (Tendulkar et al. 2014). For reference, the longest flare in the XMM-Newton observations last for ∼45 minutes and the longest optical flare seen in the MDM data (from 2013 December 28) is ∼3.6 hr long. It is unclear if this long flare is an exceptionally long and bright flare episode or that it is yet another flux mode.

As shown in Coti Zelati et al. (2014) and confirmed by our analysis, the optical emission exhibits orbital modulations, which arises due to heating of the face of the secondary star. Sporadic intense flares correlated with the X-ray ones are also seen in both the UV and optical. Neither the UVW1 or B filter photometric light curves exhibit large-amplitude, rapid flux mode switching seen in X-rays.

Based on the multi-wavelength data we have presented here, we can investigate the emission properties of PSR J1023+0038 in its present accreting state spanning the electromagnetic spectrum from the radio to the GeV γ-ray range. Figure 15 shows the spectrum of PSR J1023+0038 across nearly 18 decades of photon frequency. The X-ray spectrum is separated into the three flux modes. In Tendulkar et al. (2014) it was established that the hard X-ray emission observed with NuSTAR exhibits the same power-law spectrum and mode switching behavior as the soft X-rays. For reference, here we show the time-averaged hard X-ray spectrum that is also dominated by the high mode, as well as the Fermi LAT flux measurement for PSR J1023+0038 in the LMXB state from Stappers et al. (2014).

Zoom In Zoom Out Reset image size

One important aspect of the observed emission is how the X-ray spectrum in the low, high, and flare modes extrapolates to the UV and optical range. From Figure 15 it is evident that the emission process responsible for the X-ray flux cannot account for the majority of the optical emission as observed with XMM-Newton OM (even if we remove the sinusoidal modulations). This provides a plausible explanation for the absence of large-amplitude variations in the B filter data comparable in amplitude to the X-ray mode changes. The extrapolated high mode spectrum in the B filter contributes ≲10% to the average optical flux. The high mode spectrum can produce ∼50% of the observed emission in the Swift UVOT UVW2 filter, while in the UVW1 filter the contribution of the high mode declines to ∼30%. The increase in UVW2 emission that appears to coincide with a low-high mode transition and the lack of comparable brightness variations in the UVW1 band during X-ray mode switches are consistent with these findings.

The extrapolated flare mode spectrum falls well below the flux of the optical flares, suggesting that the X-ray and optical flares are possibly not generated by the same process. Energetically, the peak optical luminosity in the brightest flare (∼1 × 10³⁴ erg s⁻¹) is well below the peak X-ray luminosity (∼3 × 10³⁴ erg s⁻¹ for 0.3–10 keV or ∼6 × 10³⁴ erg s⁻¹ for 0.3–79 keV if we take into account the NuSTAR data) so it is plausible that reprocessing at a location other than the site of X-ray emission produces the optical flares. However, reprocessing cannot easily account for the puzzling instances where the optical peak precedes the X-ray one. The projected size of the binary is ∼4–5 light s so these significant time offsets cannot be attributed to light travel time differences due to emission from different regions of the system. A possible explanation may involve a complex emission spectrum that peaks in both the optical and X-rays.

11.1. Radio Pulsar Enshrouding or Quenching?

11.2. Pulsar Mode Switching?

11.3. Comparison with Other Systems

11.4. Previous Interpretations

11.5. The Accretion Physics of PSR J1023+0038

We have presented extensive multi-wavelength observations of PSR J1023+0038 in its LMXB state. The wealth of data offer unique insight into the behavior of accreting MSPs in the quiescent regime. The system appears to exhibit short-lived but very frequent and exceptionally stable episodes of chanelled accretion onto the stellar magnetic poles. This activity appears to be limited to a remarkably narrow luminosity range around ∼3 × 10³³ erg s⁻¹. The frequent and rapid switches to a second, lower luminosity mode with ∼5 × 10³² erg s⁻¹ during which no pulsations are observed, can be interpreted as being non-accreting intervals. At least a portion of the occasional X-ray/optical flares may be produced by enhanced, spasmodic accretion due to some form of disk instability or deformation of the magnetosphere. Although certain propeller and trapped disk models predict qualitatively similar behavior, none are fully consistent with the observed properties of PSR J1023+0038, especially the accretion-induced pulsations produced at a luminosity ∼100 times lower than previosuly observed in an AMXP.

The episodic reappearance of an accretion disk in PSR J1023+0038 implies that the current disk will eventually recede and the binary will revert to its radio pulsar state yet again. In such an event, it is important to identify any differences in the system between the pre and post LMXB state periods.

Further insight into this peculiar system can be gained with additional observations, especially contemporaneous X-ray and continuum radio observations over a long time span, to conclusively establish any correlated variability between the two bands. Perhaps most importantly, it is crucial to monitor the long-term spin behavior of PSR J1023+0038 in order to establish what kinds of torques are being imparted onto the neutron star, which may have profound implications for understanding of accretion onto highly magnetized objects, in general.

A.M.A. and J.W.T.H. acknowledge support from a Vrije Competitie grant from NWO. A.T.D. acknowledges support from an NWO Veni Fellowship. J.W.T.H. and A.P. acknowledge support from NWO Vidi grants. J.W.T.H. also acknowledges funding from an ERC Starting Grant "DRAGNET" (337062). A portion of the results presented was based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. This work was based in part on observations obtained at the MDM Observatory, operated by Dartmouth College, Columbia University, Ohio State University, Ohio University, and the University of Michigan. This research is based in part on observations made with ESO Telescopes at the Paranal Observatory under programme ID 292-5011. The WSRT is operated by ASTRON (Netherlands Institute for Radio Astronomy) with support from The Netherlands Foundation for Scientific Research. Access to the Lovell Telescope is supported through an STFC consolidated grant. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The EVN (http://www.evlbi.org) is a joint facility of European, Chinese, South African, and other radio astronomy institutes funded by their national research councils. LOFAR, the Low Frequency Array designed and constructed by ASTRON, has facilities in several countries, that are owned by various parties (each with their own funding sources), and that are collectively operated by the International LOFAR Telescope (ILT) foundation under a joint scientific policy. This research has made use of the NASA Astrophysics Data System (ADS).

Facilities: XMM, Swift