Abstract
PSR J1023+0038 is the first millisecond pulsar discovered to pulsate in the visible band; such a detection took place when the pulsar was surrounded by an accretion disk and also showed X-ray pulsations. We report on the first high time resolution observational campaign of this transitional pulsar in the disk state, using simultaneous observations in the optical (Telescopio Nazionale Galileo, Nordic Optical Telescope, Telescopi Joan Oró), X-ray (XMM-Newton, NuSTAR, NICER), infrared (Gran Telescopio Canarias), and UV (Swift) bands. Optical and X-ray pulsations were detected simultaneously in the X-ray high-intensity mode in which the source spends ∼70% of the time, and both disappeared in the low mode, indicating a common underlying physical mechanism. In addition, optical and X-ray pulses were emitted within a few kilometers and had similar pulse shapes and distributions of the pulsed flux density compatible with a power-law relation Fν ∝ ν−0.7 connecting the optical and the 0.3–45 keV X-ray band. Optical pulses were also detected during flares with a pulsed flux reduced by one-third with respect to the high mode; the lack of a simultaneous detection of X-ray pulses is compatible with the lower photon statistics. We show that magnetically channeled accretion of plasma onto the surface of the neutron star cannot account for the optical pulsed luminosity (∼1031 erg s−1). On the other hand, magnetospheric rotation-powered pulsar emission would require an extremely efficient conversion of spin-down power into pulsed optical and X-ray emission. We then propose that optical and X-ray pulses are instead produced by synchrotron emission from the intrabinary shock that forms where a striped pulsar wind meets the accretion disk, within a few light cylinder radii away, ∼100 km, from the pulsar.
1. Introduction
Transitional millisecond pulsars are fast spinning, weakly magnetic (B* ≈ 108 G) neutron stars (NSs) with a low-mass () companion star, that undergo transitions between distinct emission regimes over a timescale of less than a couple of weeks. During bright X-ray outbursts (LX ≳ 1036 erg s−1), they behave like accreting millisecond pulsars (Wijnands & van der Klis 1998; see Patruno & Watts 2012; Campana & Di Salvo 2018 for reviews), accreting matter transferred by the donor through a disk and emitting X-ray pulsations, due to the channeling of the plasma inflow onto the magnetic poles. When accretion stops (LX ≲ 1032 erg s−1), they behave as redback pulsars (D'Amico et al. 2001; Roberts et al. 2018); the rotation of the NS magnetic field powers radio and high-energy (X-rays, gamma-rays) pulsed emission, as well as a relativistic wind that shocks off the matter transferred by the companion close to the inner Lagrangian point of the binary and ejects the matter from the system (see, e.g., Burderi et al. 2001). IGR J18245–2452/PSR M28-I underwent a clear transition between these two regimes in 2013 (Papitto et al. 2013; Ferrigno et al. 2014).
State transitions from two more millisecond pulsars have been observed so far—PSR J1023+0038 (Archibald et al. 2009; Patruno et al. 2014; Stappers et al. 2014) and XSS J12270–4859 (de Martino et al. 2010, 2013; Bassa et al. 2014). However, the accretion disk state of these sources was peculiar; it lasted almost a decade and had an X-ray luminosity much lower ( erg s−1) than that usually shown by low-mass X-ray binaries (Linares 2014). The X-ray light curve repeatedly showed transitions on a timescale of ∼10 s between two intensity modes characterized by a definite value of the luminosity; these were dubbed high (LX ≃ 7 × 1033 erg s−1) and low (LX ≃ 1033 erg s−1) modes (Bogdanov et al. 2015; Campana et al. 2016). X-ray flares reaching up to a few times 1034 erg s−1 were also observed. Coherent X-ray pulsations were detected only during the high mode (Archibald et al. 2015; Papitto et al. 2015) and interpreted in terms of magnetically channeled accretion onto the magnetic poles, even though the low X-ray luminosity at which they were observed should not allow matter to overcome the centrifugal barrier due to pulsar rotation (but see also Bozzo et al. 2018). Moreover, the spin-down rate of PSR J1023+0038 during the disk state (Jaodand et al. 2016, J16 in the following) was close to the value taken during the rotation-powered radio pulsar state, and its modulus is lower by at least one order of magnitude than that expected if accretion or propeller ejection of matter takes place. A relatively bright continuous radio emission with a flat spectral shape was detected and interpreted as a compact self-absorbed synchrotron jet (Deller et al. 2015), with radio flares occurring during the X-ray low mode indicating ejection of optically thin plasmoid (Bogdanov et al. 2018). Flares and flickering reminiscent of the high-/low-mode transitions were also seen in the optical (Shahbaz et al. 2015; Kennedy et al. 2018; Papitto et al. 2018) and in the infrared (Hakala & Kajava 2018; Shahbaz et al. 2018) bands; the optical light was polarized at a ∼1% degree (Baglio et al. 2016), to a different extent during the various modes (Hakala & Kajava 2018). The appearance of a disk in transitional millisecond pulsars in such a subluminous state was also accompanied by a factor of a few increase in the GeV gamma-ray luminosity (Torres et al. 2017), while only upper limits were placed on the TeV flux (Aliu et al. 2016).
This complex phenomenology led to a flurry of different interpretations either based on the emission of a rotation-powered pulsar enshrouded by disk matter (Coti Zelati et al. 2014; Li et al. 2014; Takata et al. 2014), a pulsar that propels away disk matter (Papitto et al. 2014; Papitto & Torres 2015), or a pulsar accreting mass at a very low rate from disk trapped near corotation (D'Angelo & Spruit 2012). The possibility that switching between the high and the low X-ray modes marked changes between a radio pulsar and a propeller regime was also proposed (Linares 2014; Campana et al. 2016; Coti Zelati et al. 2018).
Recently, Ambrosino et al. (2017) discovered optical pulsations at the spin period of PSR J1023+0038, produced by a region a few tens of kilometers away from the NS. This finding was interpreted by the authors as a strong indication that a rotation-powered pulsar was active in the system, even when the accretion disk was present. In fact, the pulsed luminosity observed in the visible band was too large to be produced by the reprocessing of accretion-powered X-ray emission or cyclotron emission by electrons in the accretion columns above the pulsar polar caps. Here we report on the first high time resolution multiwavelength observational campaign of PSR J1023+0038 in the disk state aimed at exploring the relation between optical and X-ray pulses, and their properties in the different intensity modes. The campaign was performed on 2017 May 23–24 and involved simultaneous high time resolution observations by the fast optical photometer SiFAP mounted at the INAF Telescopio Nazionale Galileo (TNG), X-ray instruments on board XMM-Newton and NuSTAR, and the Canarias InfraRed Camera Experiment (CIRCE) at the Gran Telescopio Canarias (GTC). Optical, UV, and X-ray spectra and images were also obtained thanks to the Nordic Optical Telescope (NOT), Telescopi Joan Oró (TJO), and Neil Gehrels Swift Observatory observations. Further high time resolution optical observations were performed by TNG/SiFAP on 2017 December 20 and were analyzed with observations performed by the X-ray NICER mission a few weeks earlier.
2. Observations
Table 1 lists the observations analyzed and discussed in this paper. We give details on the analysis of the different data sets in the following.
Table 1. Log of the Observations of PSR J1023+0038
Telescope/Instrument | MJD Start Timea | Exposure (s) | Band | Mode/Magnitude |
---|---|---|---|---|
2017 May 23 | ||||
NuSTAR/FPMA-FPMB | 57896.035995 | 82514.0 | 3–79 keV | ⋯ |
Nordic Optical Telescope/ALFOSC | 57896.896694 | 3600.0 | 440–695 nm | grism#19 |
XMM-Newton/MOS1 (XMM1) | 57896.905272 | 24651.7 | 0.3–10 keV | Small window |
XMM-Newton/MOS2 (XMM1) | 57896.905765 | 24613.9 | 0.3–10 keV | Small window |
Telescopi Joan Oró/MEIA2 | 57896.907924 | 300.0 | Johnson V | 16.85 ± 0.01 |
Swift/UVOT | 57896.915625 | 1710.7 | UVW1 | Imaging+Event |
Swift/XRT | 57896.915579 | 1721.6 | 0.3–10 keV | Photon counting |
XMM-Newton/EPIC-pn (XMM1) | 57896.929398 | 24914.0 | 0.3–10 keV | Timing |
GTC/CIRCE | 57896.930555 | 4800.0 | Ks | Fast imaging |
TNG/SiFAP (TNG1) | 57896.970058 | 3297.7 | white filter | Fast timing |
2017 May 24 | ||||
XMM-Newton/EPIC-pn (XMM2) | 57897.739274 | 23413.0 | 0.3–10 keV | Timing |
XMM-Newton/MOS1 (XMM2) | 57897.715133 | 23150.6 | 0.3–10 keV | Small window |
XMM-Newton/MOS2 (XMM2) | 57897.715642 | 23112.8 | 0.3–10 keV | Small window |
TNG/SiFAP (TNG2) | 57897.890802 | 8397.2 | white filter | Fast timing |
Telescopi Joan Oró/MEIA2 | 57897.896660 | 600.0 | Johnson U | 18.1 ± 0.1 |
Telescopi Joan Oró/MEIA2 | 57897.903954 | 200.0 | Johnson B | 17.10 ± 0.03 |
Telescopi Joan Oró/MEIA2 | 57897.906622 | 200.0 | Johnson V | 16.64 ± 0.02 |
Telescopi Joan Oró/MEIA2 | 57897.909287 | 200.0 | Cousins R | 16.31 ± 0.02 |
Telescopi Joan Oró/MEIA2 | 57897.911955 | 200.0 | Cousins I | 15.91 ± 0.02 |
2017 Dec 2–20 | ||||
NICER (1034060101) | 58089.019028 | 3070.0 | 0.2–12 keV | ⋯ |
NICER (1034060102) | 58090.044676 | 1963.0 | 0.2–12 keV | ⋯ |
NICER (1034060103) | 58093.393750 | 1144.0 | 0.2–12 keV | ⋯ |
TNG/SiFAP (TNG3-B) | 58107.126433 | 3298.3 | Johnson B | Fast timing |
TNG/SiFAP (TNG3-V) | 58107.179563 | 3298.3 | Johnson V | Fast timing |
TNG/SiFAP(TNG3-R) | 58107.235471 | 3298.3 | Johnson R | Fast timing |
Note.
aBarycentric Dynamical Time at exposure start.Download table as: ASCIITypeset image
2.1. X-Ray Observations
2.1.1. XMM-Newton/EPIC
We analyzed XMM-Newton Target of Opportunity (ToO) observations of PSR J1023+0038 performed on 2017 May 23 (ID: 0794580801, XMM1 in the following) and May 24 (ID: 0794580901, XMM2 in the following) in the Discretionary Time of the Project Scientist (PIs: Papitto, Stella). We used SAS (Science Analysis Software) v.16.1.0 to reduce the data. We transformed the photon arrival times observed by XMM-Newton as if they were observed at the line of nodes of the solar system barycenter using the source position derived by Deller et al. (2012), R.A. = 10:23:47.687198(2), decl. = 00:38:40.84551(4); the JPL ephemerides DE405; and the barycen tool. We used the same parameters to correct arrival times observed by other instruments considered in this paper. We discarded the first 3.5 ks of XMM1 data from the analysis because high flaring particle background—apparent from the 10–12 keV light curve—contaminated data and prevented identification of X-ray modes. In both observations, the EPIC-pn was operated with a time resolution of 29.5 μs (timing mode) and a thin optical blocking filter. We defined source and background regions with coordinates RAWX = 27–47 and RAWY = 3–5, respectively, and retained good events characterized by a single or a double pattern. EPIC MOS1 and MOS2 cameras observed the target in small window mode with a time resolution of 0.3 s and a thin optical blocking filter. We extracted source photons falling within a circular region centered on the source position with a 35'' radius, and background photons from a 100'' wide, source-free circular region of one of the outer CCDs; we retained good events with patterns as complex as quadruples. We created background-subtracted light curves with the task epiclccorr. Redistribution matrices and ancillary response files were computed using rmfgen and arfgen, and the spectra rebinned to have at least 25 counts per bin and no more than three bins per resolution element in the 0.3–10 keV energy band.
2.2. NuSTAR
We analyzed the NuSTAR ToO observation of PSR J1023+0038 performed on 2017 May 23 (ID: 80201028002, starting at 00:26:09 UTC and lasting 160 ks, for a total exposure time of 82.5 ks; PI: Papitto). We reduced the observation by performing standard screening and filtering of the events with NuSTAR data analysis package (NUSTARDAS) version v.1.8.0 with CALDB 20170126. We selected source and background events from circular regions of 55'' radius centered at the source location and in a source-free region away from the source, respectively.
2.3. NICER
We present the analysis of NICER observations of PSR J1023+0038 performed on 2017 December 2 (ID: 1034060101), December 3 (ID: 1034060102), and December 6 (ID: 1034060103). The events across the 0.2–12 keV band (Gendreau et al. 2012) were processed and screened using HEASOFT version 6.24 and NICERDAS version 4.0.
2.4. Swift/X-Ray Telescope (XRT)
We consider the Swift observations of PSR J1023+0038 that started on 2017 May 23 at 21:53 (UT; ID: 00033012149) with an exposure of 1.7 ks. We reduced data obtained with the XRT in photon-counting mode using the HEASoft tool xrtpipeline, extracted light curves and spectra with xselect from a circle with a 47'' radius centered on the source position, and produced ancillary response files using xrtmkarf. The XRT observed a variable count rate between 0.05 and 0.4 s−1. The 0.3–10 keV spectrum could be described by an absorbed power law with absorption column fixed to NH = 5 × 1020 cm−2 and photon index Γ = 1.6 ± 0.2, giving an unabsorbed 0.3–10 keV flux of (1.0 ± 0.1) × 10−11 erg cm−2 s−1, typical for the disk state of PSR J1023+0038.
2.5. Optical/UV Observations
We give details on the analysis of the different optical/UV data sets in the following.
2.6. TNG
We observed PSR J1023+0038 with the SiFAP fast optical photometer (Meddi et al. 2012; Ambrosino et al. 2016, 2017) mounted at the 3.6 m TNG starting on 2017 May 23 at 23:17 (TNG1, overlapping for 3.0 ks with XMM1), on 2017 May 24 at 21:21 (TNG2, overlapping for 8.0 ks with XMM2), and on 2017 December 20 at 03:02 (TNG3). The three observations were performed during the Director's Discretionary Time. We observed PSR J1023+0038 using the channel of SiFAP that ensured the maximum possible time resolution (25 ns), and a V = 15.6 mag reference star UCAC4 454–048424 (Zacharias et al. 2013) with the second channel which operated with a time resolution of 1 ms during TNG1, and 5 ms during TNG2 and TNG3. The size of the on-source channel of SiFAP is 0.13 × 0.13 cm2, which at the F/11 focus of the TNG corresponds to nearly 7 × 7 arcsec2. In TNG1 and TNG2, we used a white filter covering the 320–900 nm band and maximum between 400 and 600 nm (i.e., roughly corresponding to the B and V Johnson filters; see Supplementary Figure 1 in Ambrosino et al. 2017). During TNG3, we observed the source with Johnson B (λeff = 445 nm, ΔλFWHM = 94 nm), V (λeff = 551 nm, ΔλFWHM = 88 nm), and R (λeff = 658 nm, ΔλFWHM = 138 nm) filters, for 3.3 ks each; filters could not be used on the reference star, due to problems with the instrumental setup. We estimated the background by tilting the pointing direction of the telescope by a few tens of arcseconds for ∼100 s; this was done four times during TNG2 (roughly every half an hour) and once during TNG1 (at the end of the exposure), and during each of the three filtered observations of TNG3. The background count rate increased during TNG2 as the elevation of the source over the horizon decreased and the contamination by diffuse light was correspondingly larger; we then evaluated the background contribution by fitting the count rate observed during the four intervals with a quadratic polynomial. During TNG1, the background was estimated only at the end of the exposure, and considering that its contribution increases as the source declines over the horizon, it is almost certainly larger than the actual value that affected the observation of the source.
The SiFAP clock showed drifts with respect to the actual time measured by two Global Positioning System (GPS) pulse-per-second (PPS) signals that were used to mark the start and stop times of each observation. For this reason, the total elapsed time by the clock exceeded the GPS time by ΔtSiFAP − ΔtGPS = 0.84 and 4.23 ms during TNG1 and TNG2, each lasting ΔtGPS = 3300.0 and 8400.0 s, respectively. During TNG3, the SiFAP clock lagged the GPS signal by 3.89, 3.13, and 3.08 ms during the exposures with the B, V, and R filters, respectively, each lasting . As no further information on the dependence of the drift on the various parameters that affect the photometer operations (e.g., temperature, count rate) is available, the best possible guess is that the drift rate was constant. Following Ambrosino et al. (2017), we corrected the arrival times recorded by SiFAP using the relation . Subsequently, we used the software tempo2 (Hobbs et al. 2006) to correct the photon arrival times to the solar system barycenter, using the position of PSR J1023+0038 reported in Section 2.1.1 and the geocentric location of the TNG (X = 5,327,447.481, Y = −1,719,594.927, Z = 3,051,174.666), along with the JPL ephemerides DE405. During the 3.3 ks exposure of TNG1, we measured an average count rate of ∼2.59 × 104, 3.70 × 104, and 1.75 × 104 s−1 from the source, the reference star, and the sky background, respectively. Values of 2.95 × 104, 4.19 × 104, and 1.11 × 104 s−1 were measured for the same quantities during TNG2.
2.7. XMM-Newton/OM
The Optical Monitor (OM) on board XMM-Newton observed the source during XMM1 and XMM2 with a B filter (λeff = 450 nm) and using the fast mode, which has a time resolution of 0.5 s. We extracted source photons from a 6 pixel wide circle (corresponding to 29) and background from an annulus with inner and outer radii of 7.2 and 9.9 pixels, respectively. Observed count rates, ROM, were converted to flux and magnitude units using the relations18
erg cm−2 s−1 Å−1 and mag = 19.266–2.5 log10 ROM.
2.8. Swift/Ultraviolet Optical Telescope (UVOT)
The UVOT on board Swift observed PSR J1023+0038 with the UVW1 filter (λc = 260 nm) for 674.1 s, starting on 2017 May 23 at 21:58. We performed aperture photometry by using a circle with radius of 5'' around the source position and extracting the background from a nearby source-free region. We measured an average flux of 16.05 ± 0.03(stat) ± 0.03(sys) mag (Vega system) from the source with the tool uvotsource, corresponding to a flux density of (1.5 ± 0.1) × 10−15 erg cm−2 s−1, or (3.38 ± 0.21) mJy.
2.9. Nordic Optical Telescope
We performed a ToO (target of opportunity) observation of PSR J1023+0038 (PI: Papitto) with the 2.5 m NOT starting on 2017 May 23 (21:28 UTC), taking four 900 s long spectra with the Andalucia Faint Object Spectrograph and Camera (ALFOSC), equipped with the 05 slit (appropriate to the actual seeing of 0
7), and grism #19 (440–695 nm). Subsequently, five photometric images were taken with the SDSS u', g', r', i', and z' filters with an exposure of 60 s (except for a 120 s exposure for the u' image).
We reduced data using standard IRAF (Tody 1986) tasks such as bias and flat-field correction, cosmic-ray cleaning, wavelength calibration, extraction of spectra from science frames using the optimized method by Horne (1986), and flux calibration. The wavelength calibrations made use of He and Ne arc lamps as a reference, while the flux calibration was made using the simultaneous observation performed by TJO in the Johnson V band (V = 16.85 ± 0.02 mag, corresponding to Fλ = (6.9 ± 0.1) × 10−16 erg cm−2 s−1 Å−1; see Section 2.10) and the La Palma standard extinction curve. The spectrum was first extracted for each differential image, then the four extracted spectra were combined together.
We used Starlink GAIA v.4.4.8 to perform aperture photometry on the five images taken by ALFOSC, using the nearby stars UCAC4 454–048421 and UCAC4 454–048418 (Zacharias et al. 2013) to calibrate the magnitude scale, and an aperture between 3'' and 5'' depending on the filter, and obtained the following magnitudes, u' = 17.43(1), g' = 16.836(4), r' = 16.625(4), i' = 16.488(5), z' = 16.370(8).
2.10. Telescopi Joan Oró
TJO is a robotic 80 cm telescope located at the Observatori Astronòmic del Montsec (Catalunya, Spain). We observed the field around PSR J1023+0038 on 2017 May 23 and 24 (see Table 1 for details) in the context of an observing program to monitor transitional millisecond pulsars (ID: p153, PI: Papitto). We used the photometric imaging camera MEIA2, equipped with Johnson U, B, and V filters and Cousins R and I filters. We created average dark and bias frames and flat-fielded images using the ESO eclipse package v.5.0–019
and performed aperture photometry with Starlink GAIA v.4.4.8 using an aperture of 12 pixels (equivalent to 43). We used the same nearby reference stars used to calibrate NOT images (see Section 2.9) and obtained magnitudes reported in the rightmost column of Table 1.
2.11. Gran Telescopio Canarias
Fast near-Infrared imaging of PSR J1023+0038 was carried out with ToO observations (PI: Rea) on the night of 2017 May 23 using CIRCE (Eikenberry et al. 2018) on the 10.4 m GTC. We configured the detector in fast imaging mode using a window size of 2048 × 1366 with the plate scale of 01 pixel−1. We obtained 475 images with 4.9 s exposure time using the Ks filter from 57,896.93 to 57,896.98 MJD. Due to fast variations in the infrared sky, we dithered the telescope every five images with a five-point dither pattern. Appropriate flat and dark frames were obtained during the twilight and at the end of the night, respectively.
Reduction of the CIRCE data was carried out using the SuperFATBOY data reduction pipeline. We applied the standard procedures in the following order: linearity correction, dark subtraction, and dividing by master flat. Infrared sky background was obtained separately for each dither pattern consisting of 25 images. In the end, we interpolated over bad pixels and cosmic-ray events, and binned the images by two pixels to improve the signal-to-noise ratio.
Extracting the photometry, we roughly aligned all of the images. Then, we further adjusted the location of each source (PSR J1023+0038 and three reference stars) with a Gaussian fit. We extracted the source counts from an aperture of 8 pixel radius and determined the sky background from an annulus with 1.5–2 times the aperture size.
3. Data Analysis
3.1. Light Curves
3.1.1. X-Ray Light Curve
We built an EPIC 0.3–10 keV light curve with a time resolution of 10 s, summing the background-subtracted light curves of the three EPIC cameras: pn, MOS1, and MOS2. We adopted the definition of X-ray modes of Bogdanov et al. (2015) and considered low-mode intervals with a count rate lower than 3.1 s−1, flaring intervals with a count rate larger than 11 s−1, and high-mode intervals with a count rate in between these thresholds. We also adopted the bistable comparator defined by Bogdanov et al. (2015) and defined an intermediate gray area ranging from 2.1 to 4.1 s−1; we did not flag this as a transition between the high mode and the low mode (or vice versa) when the count rate varied from the high-mode region to the gray area and then returned back to the high-mode region, but considered the whole interval as high mode. Figure 1 shows the light curves observed during XMM1 (panel (a)) and XMM2 (panel (c)), respectively. Low-, high-, and flaring mode intervals are plotted in red, blue, and green points, respectively. Panels (b) and (d) of Figure 1 show the OM optical light curve during XMM1 and XMM2, respectively.
Figure 1. 0.3–10 keV EPIC light curve observed during XMM1 (panel (a)) and XMM2 (panel (c)), binned every 10 s. Blue, red and green points indicate high-, low-, and flaring mode intervals (see the text for their definition). Panels (b) and (d) show the light curve observed by the Optical Monitor on board XMM-Newton and binned every 30 s. Horizontal bars indicate the intervals of simultaneous observations performed by other instruments.
Download figure:
Standard image High-resolution imageFigure 2 shows the sum of the 3–79 keV light curves observed by two NuSTAR modules, binned at a time resolution of 100 s. We identified the transition using a light curve sampled in 10 s long time bins, and setting the threshold between the low (red points) and the high (blue points) mode at a count rate of 0.15 s−1, whereas the flaring (green points) took place at a rate higher than 0.9 s−1. These thresholds are similar to those determined by Tendulkar et al. (2014) and Coti Zelati et al. (2018); we set them by requiring that the modes determined from the NuSTAR light curve would be the same as those observed by XMM-Newton during the 11.1 ks long overlap occurring with XMM1. Due to the different photons statistics, an exact correspondence could not be found; the thresholds were conservatively set in order to avoid the NuSTAR intervals deemed as flaring or low being contaminated by high-mode emission.
Figure 2. 3–79 keV light curve observed by NuSTAR binned at 100 s intervals. Red, blue, and green points mark the low, high, and flaring modes, respectively (see the text for their definition). The horizontal bar indicates the interval covered by the XMM1 observation.
Download figure:
Standard image High-resolution image3.1.2. Correlated Optical/X-Ray Variability
The top panel of Figure 3 shows the differential TNG/SiFAP optical light curve defined as the ratio between the background-subtracted count rate observed from PSR J1023+0038 and the reference star employed at the TNG, together with the simultaneous XMM-Newton/EPIC X-ray light curves in the respective bottom panels. Both light curves were binned every 10 s. Correlation between flares is evident, while there is no clear optical analog of the high-/low-mode transitions observed in the X-ray band.
Figure 3. Differential TNG/SiFAP optical light curve of PSR J1023+0038 (top panel) and simultaneous XMM-Newton/EPIC coverage (bottom panel). Blue, red, and green points indicate high-, low-, and flaring mode intervals (see text for their definition).
Download figure:
Standard image High-resolution imageIn order to explore the degree of correlation between the X-ray and optical variability, we considered the light curves observed by XMM-Newton/EPIC and TNG/SiFAP binned with a shorter timescale. The three intensity modes defined in Section 3.1.1 introduce a strong nonstationarity in the light curves. The cross-correlation function (CCF) requires stationarity; therefore, we calculated a CCF for each mode. The intervals were selected based on the X-ray behavior and were chosen carefully so as to exclude transitions. The CCFs were calculated using the HEASoft tool crosscor with a time resolution ranging from 1 to 2.5 s over 64 s intervals of the flaring, high, and low modes; they are plotted in Figure 4 with magenta, green, and blue symbols, respectively. The inset of Figure 4 shows the CCFs evaluated on a shorter timescale, ranging from 25 to 100 ms.
Figure 4. CCFs between the optical (TNG/SiFAP) and X-ray (XMM-Newton/EPIC) light curve observed during the flaring (magenta points), high (green points), and low (blue points) modes. A time resolution ranging between 1 and 2.5 s was used in the different modes, and 64 s long intervals were averaged. The inset shows the CCF evaluated on a shorter timescale (ranging between 25 and 100 ms).
Download figure:
Standard image High-resolution imageOur analysis shows that X-ray and optical emissions were clearly correlated in each of the three modes, with a somewhat similar behavior and degree of correlation. Namely, on timescales longer than a second (see Figure 4), the optical variability showed a range of delays with respect to the X-ray variability, with a reflection-like CCF shape (O'Brien et al. 2002) that rises sharply in 2–3 s, peaks around zero, and decays slowly toward positive delays (i.e., optical lagging X-rays), reaching ∼10–20 s. On shorter timescales (see the inset of Figure 4), the optical variability appeared instead to be correlated with the X-ray variability in the flaring mode with no evidence for delays. The CCF hints to a correlation also in the high mode, with no evidence for delays, while in the low mode the analysis is hampered by the lower statistics.
We note, however, that the shape of the CCFs depends strongly on the parameters chosen to calculate them, i.e., on the time resolution and on the length of the intervals over which the fast Fourier transforms were calculated. More importantly, when calculating the CCFs in the time domain (which by definition implies the use of nonstrictly simultaneous data intervals), we obtain different results. Finally, notwithstanding the rather low statistics, we have marginal evidence for the CCFs being variable in time. If confirmed, all this would imply that the variability is nonstationary even within a (X-ray) defined mode.
3.1.3. Correlated Infrared/X-Ray Variability
The infrared light curve accumulated by CIRCE at GTC overlaps for 4.3 ks with the EPIC-pn exposure in XMM1. Figure 5 shows the light curves observed in the infrared and X-ray bands with magenta and green symbols, respectively. Note that high particle background flaring affects the EPIC-pn light curve before MJD 57,896.967, i.e., during the first 2.7 ks of the overlap with the CIRCE light curve, making it extremely noisy. The combination of this high background interval affecting the X-ray light curve and the short duration and low statistics of the infrared light curve prevented us from performing a quantitative study of the correlation between the X-ray and the infrared variability. However, a visual inspection of the two simultaneous curves provides clear evidence that they are correlated, with a slightly lower infrared flux during the X-ray low mode, when compared to the high mode. This is evident especially for the low mode occurring close to MJD 57,896.976 (see Figure 5).
Figure 5. Light curves observed simultaneously by CIRCE/GTC in the Ks infrared band (magenta points; the observed rates were scaled by a factor 100) and by the XMM-Newton/EPIC in the 0.3–10 keV band (green points). A time resolution of 5.88 s was used.
Download figure:
Standard image High-resolution imageAdditionally, and perhaps more interestingly, there is also possible evidence for an increase of the infrared flux right after the transition from the low to the high X-ray mode. In other words, the start of each X-ray high-mode interval might be accompanied by a modest infrared flare. However, this evidence is based on the presence of only three such flares (i.e., those at about 0.94, 0.953, and 0.978 days since MJD 57,896.0, with the latter also seen in the TNG optical light curve; see the light curves in the left panel of Figure 3 corresponding to −2 ks in the units used there,) while in two additional occasions (0.962 and 0.971), no infrared flare seems to match the start of an X-ray high-mode interval. Further and longer simultaneous observations are needed to confirm this result.
3.2. Timing Analysis
3.2.1. 2017 May Campaign: The X-Ray Pulse
In order to perform a search for pulsations in the X-ray data sets, we first corrected for the shifts of the photon arrival time caused by the orbital motion of the pulsar in the binary system. J16 showed that the epoch of the passage of the pulsar at the ascending node of the orbit T* derived from the X-ray pulsations of PSR J1023+0038 shifts by up to tens of seconds with respect to any plausible solution based on a constant orbital period derivative. We then had to determine an orbital solution valid across the time span of the observations discussed here. We corrected the two XMM-Newton/EPIC-pn time series using the semimajor axis ( lt-s) and orbital period (Porb = 17,115.5216592 s) of the J16 timing solution, and a grid of values of T* spaced by 0.125 s around the expected value. We carried out an epoch-folding search on each of the time series by sampling each period with 16 phase bins, and estimated the best epoch of passage at the ascending node (T* = 57,896.82926(1) MJD) by fitting with a Gaussian the values of the maximum pulse profile variance found in each of the periodograms we obtained. Figure 6 shows the difference ΔT* between the epoch predicted by the radio timing solution of Archibald et al. (2013; see also Table 2 in J16) and the values measured from X-ray pulsations as a function of the number of orbital cycles performed since = 54,905.97140075 MJD. We took from Table 1 of J16 the values of ΔT* measured before ≈13,200 cycles had elapsed since (red points in Figure 6), and added the last two measurements based on the analysis presented here (blue points in Figure 6; see below and Section 3.2.3). The epoch we measured of the passage of the pulsar at the ascending node in 2017 May anticipated by = 26.8(8) s the epoch predicted by the radio timing solution, confirming the ∼20–30 s shift that occurred after the 2013 June state change already reported in J16, as well as the inability to describe with a polynomial of low order (or a sinusoidal variation; see the dashed line in Figure 6) the evolution of Tasc over the years.
Figure 6. Difference between the epoch of the passage of the pulsar at the ascending node of its orbit measured from X-ray pulsations and those expected according to the radio pulsar timing solution measured by Archibald et al. (2013) and J16 as a function of the number of orbital cycles elapsed since . Red symbols refer to epochs measured by J16, blue symbols to measures from this work. Uncertainties on each measure range from 0.05 to 0.3 s. The dashed line shows a sinusoid with period equal to ≈5300 orbital cycles, unable to explain the long-term trend.
Download figure:
Standard image High-resolution imageIn previous studies (Archibald et al. 2015, J16), X-ray pulsations were observed only during the high mode with an rms amplitude of ≃8%, and no detection was obtained during the low and flaring modes, with an upper limit of 2.4% and 1.4%, respectively. We performed a pulsation search on the mode-resolved EPIC-pn time series and obtained comparable results. Blue points in Figure 7 show the epoch-folding search periodograms obtained by folding with n = 16 phase bins the high-mode intervals observed during XMM1 (top panel) and XMM2 (bottom panel). We modeled the folded pulse profiles using two sinusoidal harmonics:
where and are the average source and background count rates, respectively, and ri and ϕi (i = 1, 2) are the rms amplitude and phase of the two harmonics employed to model the pulse profile. We estimated the average period during the two observations by modeling the phases of the first and second harmonics of the pulse profile computed over 1.2 ks long intervals using the phase residual formula (see, e.g., Equation (1) of Papitto et al. 2011), in which we let only the period of the pulsations and the epoch of passage of the pulsar at the ascending node free to vary. The second harmonic turned out to have a larger amplitude and provided the most accurate measurements, PX = 1.6879874456(5) ms, compatible with the period expected according to the J16 solution, Pref = 1.687987446019(3) ms, and T* = 57,896.829262(1) MJD. Blue points in Figure 8 show the pulse profiles obtained by folding at Pref the high-mode intervals of the two EPIC-pn observations. Table 2 lists the background-subtracted rms amplitudes ri and phases ϕi (i = 1, 2) of the pulse profile; the total rms amplitude ; the 0.3–10 keV flux Fabs, which is not corrected for interstellar absorption and evaluated by modeling the observed spectrum with an absorbed power law fixing the equivalent hydrogen column density to the value measured by Coti Zelati et al. (2014; NH = 5 × 1020 cm−2); the isotropic unabsorbed luminosity L evaluated for a distance of 1.37 kpc (Deller et al. 2012); and the pulsed luminosity Lpulsed = R × L. Upper limits on the pulse amplitude are computed at the 95% confidence level. We also measured the high-mode rms amplitude in four energy bands and used these values to evaluate the pulsed flux. Cyan and blue points in Figure 13 show the total and pulsed flux in the high mode, respectively, expressed in νFν units; a dashed line indicates a νFνpulsed ∼ ν0.3 dependence, similar to that observed for the total flux.
Figure 7. Epoch-folding search periodograms of the XMM-Newton/EPIC-pn (blue points, left scale) and TNG/SiFAP (red points, right scale) time series in the high mode (HM), and TNG/SiFAP (green points, left scale) in the flaring mode (FM), observed during 2017 May 23 (top panel) and May 24 (bottom panel). Note that the scale used for the periodogram of the TNG2 observation is reduced.
Download figure:
Standard image High-resolution imageFigure 8. Pulse profiles obtained by folding the high-mode intervals in the XMM-Newton/EPIC-pn data (blue points, left y-axis scale), the high-mode intervals in the TNG/SiFAP data (red points, right y-axis scale), and the flaring mode intervals in the TNG/SiFAP data (green points, right y-axis scale, shifted arbitrarily), around Pref, with respect to the reference epoch MJD 57,896.0.
Download figure:
Standard image High-resolution imageTable 2. Properties of the X-Ray and Optical Pulses
Instrument | Band | r1 (%) | r2 (%) | ϕ1 | ϕ2 | R (%) | Fabs | L | Lpulsed |
---|---|---|---|---|---|---|---|---|---|
High mode | |||||||||
XMM-Newton/EPIC-pn | 0.3–10 keV | 2.9(3) | 8.0(3) | 0.33(2) | 0.434(3) | 8.5(4) | 13.3(1) | 30.5(2) | 2.6(1) |
NuSTARa | 3–79 keV | <1.6 | 8.4 ± 1.3 | ⋯ | ⋯ | 8.4 ± 1.3 | 22(1) | 49(2) | 4.1(7) |
TNG/SiFAPb | 320–900 nm | 0.63(2) | 0.74(2) | 0.473(7) | 0.548(3) | 0.97(2) | 4.2(2) | 11.5(4) | 0.111(5) |
Low mode | |||||||||
XMM-Newton/EPIC-pn | 0.3–10 keV | ⋯ | ⋯ | ⋯ | ⋯ | <3.0 | 2.6(1) | 6.1(2) | <0.18 |
NuSTARa | 3–79 keV | ⋯ | ⋯ | ⋯ | ⋯ | <14.0 | 2.4(4) | 5.4(9) | <0.75 |
TNG/SiFAP | 320–900 nm | <0.03 | <0.02 | ⋯ | ⋯ | <0.034 | 4.4(2) | 12.1(4) | <0.004 |
Flaring mode | |||||||||
XMM-Newton/EPIC-pn | 0.3–10 keV | ⋯ | ⋯ | ⋯ | ⋯ | <1.3 | 37.6(5) | 87.1(9) | <1.1 |
NuSTARa | 3–79 keV | ⋯ | ⋯ | ⋯ | ⋯ | <6.9 | 61(8) | 137 ± 18 | <9.5 |
TNG/SiFAP | 320–900 nm | 0.12(2) | 0.11(2) | 0.37(2) | 0.58(1) | 0.16(2) | 8.5(4) | 23(1) | 0.037(5) |
Notes. Pulse profiles were obtained by folding the respective time series in 64 phase bins (reduced to 16 in case of a low significance detection or to derive upper limits) around Pref = 1.687987446019 ms and setting Tref = 57,896.0 MJD as the reference epoch. ri and ϕi are the background-subtracted rms amplitude and phase (in cycles), respectively, of the first (i = 1) and second (i = 2) harmonic used to model the pulse profiles; R = (r12 + )1/2 is the total rms amplitude; Fabs is the flux in units of 10−12 erg cm−2 s−1 not corrected for interstellar absorption and observed in the bands listed in the second column; L is the isotropic luminosity in units of 1032 erg s−1 corrected for interstellar absorption and evaluated for a distance of 1.37 kpc (Deller et al. 2012); and Lpulsed = R × L is the pulsed luminosity evaluated using the same parameters.
aNuSTAR pulse profiles were folded around PNuSTAR = 1.687987299(1) ms, which differs significantly from Pref, due to NuSTAR clock drifts. Phases are not reported as the NuSTAR clock drifts prevent a meaningful comparison with those measured with the other instruments to be drawn. bIn TNG2, we considered only the high-mode intervals before the onset of the long flaring event, which started ∼82 ks after May 24 at 00:00 (see Figure 3 and text for details).Download table as: ASCIITypeset image
We searched the 3–80 keV NuSTAR time series for a coherent signal, computing power density spectra over 3.3 ks intervals. The average of the 28 power spectra so obtained shows an excess centered at 1184.8430(1) Hz interpreted as the second harmonic of the signal at the spin period of PSR J1023+0038. The peak in the power density spectrum is broad with a width compatible with the spurious derivative of 10−10 Hz s−1 introduced by the NuSTAR clock drift (Madsen et al. 2015) and corresponds to an rms pulse amplitude of (6.8 ± 1.3)%. To determine a more precise value of the pulsation period, we performed an epoch-folding search of the entire NuSTAR data around the corresponding fundamental frequency obtained from the analysis of the power density spectrum with steps of s, for a total of 10,001 steps. The pulse profile with the largest signal-to-noise ratio corresponds to the period PNuSTAR = 1.687987299(1) ms that differs by 1.5 × 10−10 s from Pref. We also searched for a signal in the time series restricted to the three modes. The signal was detected in the high mode with an overall rms amplitude of (8.4 ± 1.3)%, roughly constant up to 45 keV. The pulsed flux spectral distribution (see green points in Figure 13, where the total flux is also plotted in light green points) observed by NuSTAR in the high mode is compatible with the power-law relation indicated by XMM-Newton data. On the other hand, upper limits of 14.0% and 6.9% (95% confidence level) were set on the pulse amplitude during the low and flaring modes, respectively.
3.2.2. 2017 May Campaign: The Optical Pulse
We used the low-, high-, and flaring-mode time intervals determined using the EPIC light curves to select the corresponding intervals during the TNG observations. Table 2 lists the properties of the optical pulse during the three modes. We detected optical pulses at a high confidence level during the high mode (red points in Figure 7 show the periodogram) with an average, background-subtracted rms amplitude of (0.68 ± 0.02)%. Note that because the background is probably overestimated during TNG1 (see Section 2.6), the actual intrinsic pulse amplitude was likely slightly larger. The top panel of Figure 9 shows the evolution of the rms amplitude over high-mode time intervals of length ranging from ∼200 to 450 s; it attains values as high as Rmax = (1.4 ± 0.1)% and is variable although with no clear correlation with the orbital phase. Toward the end of TNG2, the pulse amplitude decreased down to a level comparable to that observed during flares (<0.2%; see below), suggesting that the entire interval from ∼82 ks after May 24 00:00 to the end of TNG2 (see right panels of Figure 3) was in the flaring mode. By considering only the high-mode intervals before the onset of such a long flaring event, the average pulse fraction observed in the high mode increased to RHS = (0.97 ± 0.02)%. Similar to the X-ray band, the optical pulse was not detected during the low mode down to an upper limit of 0.02% (95% confidence level), i.e., roughly 30 times smaller than the amplitude observed during the high mode. The flaring mode took place only during TNG2, and the optical pulse was detected at a significance of 8.3σ and with amplitude R2,flare = (0.16 ± 0.02)%, i.e., more than five times smaller than during the high mode. Considering that the net count rate observed by SiFAP/TNG during flares is about twice that in the high mode, the amplitude decrease is larger than what would be expected if the flares were a simple addition of unpulsed flux to the high-mode level. We checked that optical pulses were detected at a significance larger than 4.5σ even when flares were identified by using a higher threshold in the EPIC X-ray light curve, 20 s−1, rather than 11 s−1, the value used by Bogdanov et al. (2015) and throughout this paper.
Figure 9. Top panel: the rms amplitude of the optical pulse along the May 23 (left) and May 24 (right) observing runs during the high mode. Lower panel: the phases of the optical pulse as measured in the May 23 (left) and May 24 (right) TNG runs (blue dots and sky-blue empty squares were used for the first and second harmonic, respectively), and the phases of the X-ray pulse as measured in the two XMM-Newton runs (green dots and red empty squares were used for the first and second harmonic, respectively).
Download figure:
Standard image High-resolution imageWe performed pulse phase timing on the optical pulse profiles computed over intervals of length spanning between 200 and 500 s in the high mode. We measured an optical period Popt = 1.687987445(1) ms, a value compatible with Pref and PX, and the epoch of passage at the ascending node = 57,896.829267(6) MJD, compatible with the X-ray estimate within δT* ∼ 0.5 s. This estimate can be used to constrain the position of the region emitting the optical pulses within km in azimuthal distance along the orbit from the region emitting the X-ray pulses, where a1 is the semimajor axis of the pulsar orbit and i is the system inclination.
Red and green points in Figure 8 show the normalized and background-subtracted optical pulse profiles observed in the high and flaring modes, respectively. The optical pulse is described by two harmonics with an amplitude ratio r2/r1 ≃ 1, smaller than that of the X-ray pulse (≃3; blue points in Figure 8)). The phase of both harmonic components of the optical pulse profile lag the corresponding components of the X-ray profile; the lags of the first and second harmonic are δϕ1 = 0.14 ± 0.01 and δϕ2 = 0.112 ± 0.004, respectively. The phase difference observed during periods of strictly simultaneous observations is compatible with these estimates, indicating that the phase lag is not due to variability of the X-ray pulse profile in intervals not overlapping with the optical observations. The bottom panel of Figure 9 shows the evolution of the phase of the two harmonic components of the optical (blue dots and cyan hollow squares for the first and second harmonic, respectively) and X-ray (green dots and red hollow squares) pulse profiles over the interval of simultaneous coverage; a phase shift of ≈0.1 is always observed, with a somewhat larger timing noise shown by the optical pulse.
We measured the average optical flux Fabs in different modes by scaling the observed background-subtracted count rate by the conversion factor determined in Appendix A.1. We estimated the pulsed luminosity by scaling these flux values for the ratio between the dereddened and the absorbed flux in the 320–900 nm energy band (kder = 1.217), and by multiplying for the total rms amplitude, we obtained the pulsed luminosity shown in the rightmost column of Table 2. The values obtained in this way are listed in Table 2. A red point in Figure 13 shows the average dereddened pulsed flux observed during the high mode (5.2(4) × 10−14 erg cm−2 s−1) in νFν units, computed assuming a constant value over the 320–900 nm band.
3.2.3. 2017 December Campaign: The Optical Pulse
We preliminarily measured the epoch of passage of the pulsar at the ascending node, correcting the time of arrival of TNG3 with the semimajor axis and orbital period of the J16 timing solution, and varying a grid of values of Tasc spaced by 0.2 s around the expected value. We corrected the event arrival times of TNG3 with the preliminary value of Tasc and folded the resulting time series around Pref. We measured the pulse phase of the first and second harmonic of the pulse and determined PTNG3 = 1.687987456(12) ms (compatible with the value expected according to the J16 solution, 1.687987446177(4) ms), and = 58,107.009511(4) MJD. Blue, green, and red points in Figure 10 mark the rms amplitude (top panel) and the phases (dots show those computed on the first harmonic, hollow squares indicate phases of the second harmonic) of the pulse profiles observed during the exposures performed with the B, V, and R filters, respectively, and measured using the J16 ephemerides and the same reference epoch used to fold 2017 May data. The phase difference of ≃0.2 cycles shown by the phases computed over the first and second harmonic with respect to the values observed during 2017 May (see Table 2 and bottom panel of Figure 9) is compatible with the 0.1 cycle phase uncertainty obtained by propagating the errors on the ephemerides given by Jaodand et al. (2016) to the epoch of TNG3. A comparison between the pulse phase measured at different epochs is further hampered by the noise affecting the measured epoch of the passage at the ascending node (see Figure 6). Note that the absence of simultaneous X-ray observations prevented us from identifying transitions between the high and low modes (flares are not evident in the light curves); because pulses are absent in the low mode, the amplitudes plotted in Figure 10 are underestimated with respect to the values measured in the high mode alone. Note that the source spends on average ∼70%–80% of the time in the high mode (J16).
Figure 10. The top panel shows the rms amplitude of the optical pulse profiles observed during TNG3-B (blue points), TNG3-V (green points), and TNG3-R (red points) by folding the time series around Pref and choosing MJD 58,107.0 as the reference epoch. The bottom panel shows the phase computed over the first (filled circles) and the second harmonic (hollow squares).
Download figure:
Standard image High-resolution image3.2.4. 2017 December Campaign: The X-Ray Pulse
We analyzed three NICER observations performed in 2017 December. The bottom panel of Figure 11 shows a sample light curve observed during the first observation, where the high, low, and flaring modes can be easily recognized.
Figure 11. Pulse profile obtained by folding the high mode intervals in the NICER data with respect to the reference epoch MJD 57,896.0, using J16 ephemerides and the epoch of passage at the ascending node determined from 2017 December optical data (see Section 3.2.3), MJD (top panel). Light curve observed by NICER during observation 1034060101 (see Table 1, bottom panel).
Download figure:
Standard image High-resolution imageAdopting the orbital ephemerides from the timing analysis of the TNG3 data, we corrected the NICER photon arrival times and searched for coherent pulsations by computing power density spectra over 2 ks long intervals. A statistically significant excess at ν = 1184.8426(5) Hz is observed in the average power density spectrum, consistent with the second harmonic of the spin frequency of the source. We applied epoch-folding search techniques to the available NICER observations around the spin period extrapolated from the coherent signal detected with the Fourier analysis. We obtained the best profile for the spin period PNICER = 1.6879874457(4) ms, compatible with the value expected according to the J16 solution. The average profile is characterized by an rms amplitude R = (5.3 ± 0.7)%. Extrapolating the count rate threshold adopted for XMM-Newton/EPIC-pn data, we investigated coherent signals in the time series restricted to the three modes. Similarly to the XMM-Newton/EPIC-pn case, we did not detect any signal in the low and flaring modes. On the other hand, we detected pulsation in the high mode with amplitude of (5.6 ± 0.9)%, with a second harmonic roughly three times times stronger than the fundamental (see top panel of Figure 11).
3.3. The Optical Spectrum
Figure 12 shows the optical spectrum observed by ALFOSC/NOT, calibrated in flux using the simultaneous measurement of the magnitude by TJO (see Section 2.10). Broad, and in a few cases double-peaked, emission lines of Hα 656.3, Hβ 486.1, He i 587.6, He i 667.8, He i 501.5, He ii 468.6, and He i 447.1 nm are most prominent. Telluric contamination produced the absorption line at ≃687 nm. These features have double-peaked profiles, which are signatures of an accretion disk viewed at moderate inclination. As the continuum is blue and smooth (see Figure 2 in Wang et al. 2009), we created a template spectral shape, f1023,NOT(λ) by extrapolating the continuum observed at NOT to the 320–900 nm band, giving a flux of FJ1023,NOT ≃ 3.8 × 10−12 erg cm−2 s−1 over this band.
Figure 12. Optical (440–695 nm) spectrum of PSR J1023+0038 taken with ALFOSC at the NOT. The most prominent emission lines are marked with a vertical segment. The spectrum was normalized to match the flux density observed simultaneously by TJO in the V band (V = 16.85 ± 0.02 mag, corresponding to (6.9 ± 0.1) × 10−15 erg cm−2 s−1 nm−1 at 550 nm. During the same interval of NOT observations, the OM on board XMM-Newton observed a count rate of 6.1(1) s−1, which translates into a flux density of 7.9 × 10−15 erg cm−2 s−1 nm−1 at 450 nm.
Download figure:
Standard image High-resolution image4. Discussion
This paper presented the first high time resolution optical/X-ray/IR/UV observational campaign of PSR J1023+0038 in the disk state. Similar to other transitional millisecond pulsars in such a state, the X-ray light curve of PSR J1023+0038 shows three intensity modes (dubbed low, high, and flaring), and coherent X-ray pulsations were detected only in the high mode; they have an rms amplitude of ≃8% and are detected up to 45 keV, as first shown in this study based on NuSTAR observations.
Our simultaneous optical/X-ray observations revealed that optical pulses are also observed in the high mode. Optical pulses have an average rms amplitude of RHS = (0.97 ± 0.02)%, corresponding to a pulsed luminosity of erg s−1 (see Table 2). The spin period and the epoch of passage at the ascending node measured from the optical pulses are consistent with the values measured in the X-ray band within (0.6 ± 1.2) × 10−12 s and (0.4 ± 0.5) s, respectively. The latter estimate indicates an azimuthal distance of a few kilometers at most between the region emitting the optical and X-ray pulses, indicating that the optical and X-ray pulses are produced in the same region. Assuming a constant flux over the 320–900 nm band, Fν ∼ counts, we estimated the average pulsed flux density of . This value is compatible with the low-frequency extrapolation of the Fν ∝ ν−0.7 trend that holds for the pulsed flux density in the X-ray band (see Figure 13), suggesting that such a component could describe the pulsed spectral energy distribution over at least four decades in energy. The optical pulse amplitude was highly variable between 0.5% and 1.5% over 500 s long intervals, corresponding to a maximum pulsed luminosity of ≃ 1.6(1) × 1031 erg s−1.
Figure 13. Total and pulsed spectral energy distribution of PSR J1023+0038 in the high mode, corrected for interstellar extinction. The NOT optical spectrum (see Section 3.3) is plotted with an orange line. The average pulsed TNG/SiFAP optical flux observed in 2017 May (see Section 3.2.2) is plotted with a red square. The flux observed by Swift/UVOT UVW1 is plotted with a yellow square. The total and pulsed X-ray fluxes observed by XMM-Newton (see Section 3.2.1) are plotted in light and dark blue points, respectively. The total and pulsed X-ray fluxes observed by NuSTAR (see Section 3.2.1) are plotted in light and dark green points, respectively. Magenta squares show the average Fermi-LAT spectrum measured by Torres et al. (2017). The dashed line is a cutoff power law normalized to match the optical and X-ray observed pulsed flux.
Download figure:
Standard image High-resolution imageSimilar to X-ray pulses, optical pulses were not detected in the low mode, with an upper limit of RLS < 0.034% (95% confidence level). This corresponds to a pulsed luminosity of 4 × 1029 erg s−1, i.e., more than 25 times smaller than the value observed in the high mode.
The simultaneous detection of pulses in the optical and X-ray bands breaks down during flares. Optical pulses were detected with an rms amplitude of RFL = (0.16 ± 0.02)%, corresponding to an average pulsed luminosity of 0.37(5) × 1031 erg s−1, i.e., almost one-third than the average value observed in the high mode. X-ray pulsations remained undetected down to an amplitude of 1.3%, corresponding to a pulsed flux 2.3 times lower than that in the high mode. If the X-ray pulsed fraction decreases during flares with respect to the high mode by the same amount seen in the optical, we expect them at an amplitude of ∼1%, i.e., slightly lower than the upper limit we derived. We conclude that the nondetection of X-ray pulses during flares may result from limited photon statistics.
Both the optical and X-ray pulses were described by the sum of two harmonic components, yielding two intensity peaks every NS rotation. The ratio of the second to the first harmonic amplitude of the optical pulse r2/r1 was close to unity, lower than the value observed in the X-ray band (≃3). The first and the second harmonic of the optical pulse lagged the X-ray pulse by δτ1 = 236 ± 17 μs and δτ2 = 189 ± 7 μs, respectively. Lags of the same order were observed over intervals of a few hundred seconds in both simultaneous optical/X-ray observations analyzed here. The absolute timing accuracy of XMM-Newton is 48 μs (Martin-Carrillo et al. 2012). On the other hand, we estimated the SiFAP absolute time accuracy as better than ∼60 μs by relying on an hour-long observation of the Crab pulsar performed at the Cassini Telescope at Loiano Observatory (see the Appendix). The resulting total systematic error affecting the phase lag between the optical and X-ray pulse is <77 μs.
Optical pulses were seen at an amplitude varying between 0.5% and 2% in the Johnson B, V, and R filters. In the absence of simultaneous observations of a reference star and its X-ray counterpart, we could not measure accurately the optical spectral energy distribution.
The main result presented in this paper is that optical and X-ray pulses closely trace the repeated transitions between the X-ray high mode, in which they are both observed and whose pulsed flux is compatible with a single power-law relation , and the low mode, in which they are not detected. This strongly indicates that optical and X-ray pulses are related to the same phenomenon, something also suggested by the similar pulse shape and small phase offset. In the following, we discuss the implication of these findings for the different scenarios proposed to explain the enigmatic nature of the disk state shown by PSR J1023+0038 and other transitional millisecond pulsars.
4.1. Accretion onto the NS Surface and the Propelling of Matter
Coherent X-ray pulsations seen in the high mode of PSR J1023+0038 were first interpreted as due to magnetically channeled accretion onto the NS hot spots (Archibald et al. 2015). This interpretation was justified by the tenfold increase in the X-ray pulsed and total flux that occurred after an accretion disk formed in the system in 2013 June, as well as the similarity of the pulsed fraction, shape, and spectrum of the X-ray pulses observed from PSR J1023+0038 to those shown by accreting millisecond pulsars. Based on a similar reasoning, Papitto et al. (2015) also interpreted in terms of accretion power the X-ray pulsations shown by the transitional millisecond pulsar XSS J12270–4859 in the disk state. However, as both authors noted, interpreting such X-ray pulses as due to channeled accretion challenged the usual accretion/propeller picture for fast rotating NSs. To show this, we assume that the disk is truncated at a radius equal to ξ times the Alfvén radius,
where ξ0.5 = (ξ/0.5), is the mass accretion rate in units of 1014 g s−1 and μ26 is the magnetic dipole moment in units of 1026 G cm−3 ( for PSR J1023+0038, with the α angle between the magnetic and the spin axes, evaluated using the spin-down rate during the radio pulsar phase measured by J16 in the relation given by Spitkovsky 2006). The observed X-ray luminosity (≃8.5 × 1033 erg s−1 in the 0.3–80 keV band; Coti Zelati et al. 2018) indicates a mass accretion rate of , where m1.4 is the NS mass in units of 1.4 M⊙, and according to Equation (2), a disk radius of rm ≃ 72 ξ0.5 km. Taking ξ0.5 = 1, this value exceeds by three times the corotation radius of PSR J1023+0038 (rc = 23.8 km), and accretion onto the NS surface should be inhibited by the quick rotation of magnetic field lines at the magnetospheric boundary.
Papitto & Torres (2015) proposed that the magnetosphere would be squeezed to ∼rc if the actual accretion rate in the disk were larger ( = 5–17) than the value indicated by the X-ray luminosity. In such a case, 99% of the mass should be ejected from the inner rim of the disk by the propelling magnetosphere in order to match the relatively low X-ray flux. Roughly half of the X-ray emission and the entire gamma-ray output would be due to synchrotron self-Compton emission by electrons accelerated in the shock formed at the disk–magnetosphere boundary (see also Papitto et al. 2014). Alternatively, the disk could have fallen in a low- trapped state (D'Angelo & Spruit 2012) with its inner rim staying close to the corotation radius regardless of how low the mass accretion rate gets. In the case of PSR J1023+0038, such a model would imply that the observed X-ray luminosity traces the actual mass accretion rate and a strong outflow would not be launched.
Assuming that X-ray pulses are due to the accretion of matter onto the magnetic poles, the optical pulses of PSR J1023+0038 could be explained by cyclotron emission by electrons inflowing in magnetized accretion columns, similar to the case of accreting magnetic white dwarfs (Masters et al. 1977; Lamb & Masters 1979). Indeed, the fundamental cyclotron energy for electrons inflowing in accretion columns permeated by a magnetic field of the order of that estimated for PSR J1023+0038 (Archibald et al. 2013; J16) is Ecyc = 1.2 eV. Assuming that β is the angle between the magnetic axis and the disk plane, matter inflowing from a disk truncated at the corotation radius rC forms hot spots on the NS surface of size:
i.e., ∼10% of the NS surface, in agreement with the results of the simulations performed by Romanova et al. (2004). The corresponding typical accretion column transverse length scale is then ℓ ≃ 5 km, and the electron density in the accretion columns for a fully ionized plasma is
where μe = 1.18 is the mean molecular weight per electron, mH is the proton mass, is the freefall velocity close to the NS surface, and ℓ5 ≡ ℓ/5 km. The resulting optical depth to cyclotron absorption of accretion columns filled by plasma of such density and permeated by a B ∼ 108 G magnetic field is as large as ≈106 in the transverse direction (Trubnikov 1958). This ensures that the emission of the first few cyclotron harmonics is self-absorbed (up to roughly 10), and that the resulting spectrum is described by the Rayleigh–Jeans section of a blackbody spectrum with temperature equal to the electron temperature, Te.
However, in Ambrosino et al. (2017), we estimated the maximum cyclotron luminosity expected in the 320–900 nm band as
where νl and νh are the boundaries of the band observed by SiFAP. This value is ≈40 times lower than the observed average optical pulsed luminosity. This discrepancy holds even when taking for kTe a value of ∼100 keV, of the order of that observed from accreting millisecond pulsars (Patruno & Watts 2012 and references therein), and likely an overestimate of the temperature attained by electrons in the accretion column of a pulsar with an accretion luminosity of <1034 erg s−1. In fact, such a high electron temperature can be reached if the pressure exerted by the radiation emitted from the hot spots balances the gravitational inflow of plasma in the accretion columns and forms a shock standing off the NS surface, where the kinetic energy of the flow is converted into thermal motion of the charges. The critical luminosity to form such a shock is Lcrit ≈ 1036 ℓ5 erg s−1 (Basko & Sunyaev 1976), more than a hundred times larger than the value observed from PSR J1023+0038. Below such a value, the ions of the infalling plasma are best slowed down by Coulomb collisions with atmospheric electrons (see, e.g., Frank et al. 2002), and a temperature of the order of the effective blackbody temperature is attained, kTe ≃ 1.9 (L/Lcrit)1/4 keV. We deduce that magnetic accretion at a rate L ∼ 10−3 Lcrit is hardly capable of producing electrons hot enough to yield a sizable cyclotron emission in the optical band (e.g., Lcyc ≃ 1027 erg s−1 for kTe ≃ 0.3 keV).
According to Equation (5), the maximum observed pulsed optical luminosity erg s−1 corresponds to an unrealistically large brightness temperature of kTb = (175 ± 10) (rem/km)−2 MeV, where rem is the radius of the emission region. Even considering rem ∼ 100 km (i.e., approximately the size of the light cylinder), the X-ray luminosity that would be produced by such a thermal component ( erg s−1) would be huge. We can then safely rule out emission from hot spots on the NS surface heated by the accretion flow as the origin of optical pulses. A similar reasoning rules out the reprocessing of the X-ray emission by the inner regions of the disk, which would necessarily produce an even cooler and fainter thermal spectrum,
4.2. Rotation-powered Pulsar
An alternative possibility is that the optical pulsations of PSR J1023+0038 originate in the activity of a rotation-powered pulsar. So far, optical pulsations have been detected from five isolated high-magnetic-field young pulsars (Mignani 2011 and references therein). Models envisage that synchrotron emission of secondary electron/positron pairs accelerated in magnetospheric gaps, reconnection events, and/or the equatorial current sheet (see, e.g., Venter et al. 2018 for a recent review) give rise to nonthermal pulsed emission at optical and X-ray energies (Pacini & Salvati 1983), whereas curvature radiation accounts for the gamma-ray emission (Romani 1996). Recently, the need to use a common description of these processes, dubbed as synchro-curvature, has become evident as both effects are relevant along the particles' trajectories in the magnetosphere (Viganò et al. 2015). These mechanisms are unlikely to work if the magnetosphere is engulfed by plasma from the disk (density of ne ≃ 5 × 1015 cm−3, see Equation (4), i.e., ≈106 times the Goldreich & Julian critical density; Goldreich & Julian 1969) as gaps in the outer magnetosphere would be readily filled20 (Shvartsman 1971). Even if electron acceleration up to a Lorentz factor γ3 ≡ γ/103 = 1 were possible, charges would be stopped down by Coulomb collisions with ions and electrons of the plasma on a typical length scale of cm (see, e.g., Equation (3.35) in Frank et al. 2002). This is much smaller than the length over which electrons radiate synchrotron X-ray and optical photons in a pulsar magnetosphere, <10−5rLC ≃ 80 cm for PSR J1023+0038 after particle injection (Torres 2018). For this reason, we assume that a rotation-powered pulsar is able to work only if the matter inflow is truncated outside the light cylinder, and no accretion onto the NS surface takes place. To meet this condition, the mass accretion rate should be lower than < (see Equation (2) for rm > rLC), and the inferred luminosity lower than 2.5 × 1032 erg s−1. In this scenario, most of the observed X-ray luminosity (≃8.5 × 1033 erg s−1) would not be due to disk accretion.
However, assuming that the optical pulses of PSR J1023+0038 originate in the magnetosphere of a rotation-powered pulsar presents a number of issues. First, a very large efficiency is required to explain the conversion of up to of the spin-down power erg s−1 (Archibald et al. 2013) in the 320–900 nm optical pulsed luminosity. Values lower by at least an order of magnitude are observed from other powerful rotation-powered pulsars, including the Crab pulsar (Percival et al. 1993) and the isolated millisecond pulsar PSR J0337+1715 (Strader et al. 2016; see Figure 3 of Ambrosino et al. 2017, which compares the optical efficiency in the B band of various types of pulsars).
Second, the pulsed X-ray luminosity in the high mode was ≃2.6 × 1032 erg s−1 (Archibald et al. 2015; see also Table 2), i.e., ∼6 × 10−3 times the spin-down power. The simultaneous detection of optical and X-ray pulses in the high mode and their disappearance in the low mode means that if the former has a magnetospheric origin, the latter likely also does. The fraction of the spin-down power converted into X-ray pulses of PSR J1023+0038 would be much larger than that of almost all rotation-powered pulsars (<10−3; Possenti et al. 2002; Vink et al. 2011; see also Lee et al. 2018 for an updated analysis that suggests an average efficiency ≃10−4). Considering that rotation-powered pulsars with sinusoidal pulse profiles generally show a pulsed fraction of ∼10% (Zavlin 2007), the discrepancy is likely even larger. More importantly, when a disk was absent and radio pulses were observed, Archibald et al. (2010) detected X-ray pulsations in the 0.25–2.5 keV band with an rms amplitude of (11 ± 2)%. Even assuming that pulsations were present also in the 2.5–10 keV band and with the same amplitude (note that in that band Archibald et al. 2010 could only place a 20% upper limit on the pulse amplitude), the pulsed X-ray luminosity would have been erg s−1, i.e., 2 × 10−4 times the spin-down power. The 25 fold increase of the pulsed flux that occurred when a disk formed in the system would be very difficult to explain assuming that the rotation-powered pulsar kept working as if it were in the radio pulsar state. One case of mode switching by an isolated rotation-powered pulsar is known (Hermsen et al. 2013; Mereghetti et al. 2016), but it appears contrived that the mode change of PSR J1023+0038, which occurred when a disk formed, was not influenced by it.
The large optical and X-ray spin-down conversion efficiency needed to produce a magnetospheric emission large enough to explain the observed pulsed flux could be indeed related to the presence of the disk. Soft disk photons could enhance the pair production in the magnetosphere, yielding a brighter pulsed radiation than that in the radio pulsar state in which the system roughly behaves as if it were isolated. However, a simultaneous fit of the gamma-ray and X-ray emission of PSR J1023+0038 with models developed for rotation-powered pulsars was troublesome. We considered the synchro-curvature model developed by Torres (2018). The model was shown to be able to describe well the X-ray and gamma-ray emission of rotation-powered pulsars in terms of a few-order parameters (such as the accelerating electric field, a measure of how uniform the distribution of particles emitting toward us is, the magnetic gradient along a field line, and a normalization). Particularly, in all cases in which both energy bands displayed pulsed emission, the spectral model built out of only the gamma-ray data is close to the detected X-ray emission spectrum already, and further common analysis of both energy regimes makes for a perfect agreement (Li et al. 2018; Torres 2018). This has proven to not be the case here: we attempted to model the gamma-ray/X-ray pulsed energy distribution observed from PSR J1023+0038 (see Figure 13), but even assuming that the gamma-ray emission comes entirely from the magnetosphere (note that gamma-ray pulses have not been detected from PSR J1023+0038 in the disk state so far), varying the relevant magnetospheric parameters over the wide range used in Torres (2018) gives an X-ray and optical pulsed output lower than observed by one and three orders of magnitude, respectively. Based on the different behaviors found for PSR J1023+0038 when compared to all other pulsars studied from the synchro-curvature model, we conclude that the magnetospheric activity of a rotation-powered pulsar that works as if it were isolated (and with most of the gamma-ray radiation pulsed) is unlikely the only source of the optical/X-ray pulses observed from PSR J1023+0038.
4.3. Pulsar Wind
We suggest an alternative interpretation of the optical and X-ray pulses shown by PSR J1023+0038 in terms of synchrotron radiation emitted from the intrabinary termination shock of the pulsar wind with the accretion disk inflow at a distance krLC, with k = 1–2, i.e., just beyond the light cylinder (see Figure 14 for a schematic diagram of the geometry we have assumed; see also Veledina et al. 2019, who have presented an interpretation based on similar assumptions). For an isotropic pulsar wind, the postshock magnetic field is (Arons & Tavani 1993)
where σ is the magnetization parameter of the wind (Kennel & Coroniti 1984), which is ≫1 close to the light cylinder as the whole pulsar wind energy is carried by the electromagnetic Poynting flux (Arons 2002), and feq is a geometric factor that defines the fraction of the sky into which the pulsar wind is emitted and is unity if the wind is isotropic. For small values of k, i.e., not far from the light cylinder, the medium is permeated by such a large magnetic field that synchrotron emission is the dominant cooling mechanism for electrons accelerated at the shock. A single population of electrons with energy spectrum NE ∼ E−2.35 and cutoff at ∼2 GeV (see, e.g., Equation (35) of Lefa et al. 2012) would produce a spectral energy distribution compatible with the shape suggested by the pulsed flux measured both in the optical and the 0.3–45 keV X-ray band, νFνpulsed ∼ ν0.3 (see the dashed line in Figure 13; Martín et al. 2012). This population could result from a Fermi process with an acceleration parameter ξ ≃ 0.01 (see Equation (20) in Papitto et al. 2014; Papitto & Torres 2015). At low energies, the synchrotron emission becomes optically thick below eV (Rybicki & Lightman 1979). Because the electron density at the innermost regions of a disk truncated at k rLC with mass accretion rate < 0.3 (see Section 4.2) is / × ≃ 5 × 1012k−3/2 cm−3, the break to optically thick emission is expected below ∼0.2 eV. This ensures that the shock region (see light gray shaded region in Figure 14) is optically thin to emission in the visible band.
Figure 14. Schematic diagram of the pulsar wind scenario we propose to explain the optical and X-ray pulses observed from PSR J1023+0038. Dashed lines represent the current sheet, which expands from the light cylinder surface at rLC as an Archimedean spiral. Its arms cross the termination shock (shaded in light gray) at S1 and S2 where particles are accelerated to relativistic energies. As long as the particle energy is quickly radiated away and the size of the shock is smaller than a few light cylinder radii, two bright synchrotron emitting spots (drawn in blue) rotate at the shock surface. An observer will see pulses of radiation because the intensity received from S1 is modulated by the angle under which the spot is seen, and the emission coming from S2 is absorbed by the optically thick disk inflow (shaded in dark gray).
Download figure:
Standard image High-resolution imageClose to the light cylinder, the spin-down energy of the pulsar is transported outward by the electromagnetic field. In the striped wind model, the magnetic field configuration has the shape of two monopoles of opposite polarity which join at the equatorial plane (Bogovalov 1999), and a current sheet forms along such a plane where the field changes polarity. In the oblique rotator case, the rotation of the pulsar introduces oscillations in the the current sheet, which expands as an Archimedean spiral with arms separated by πrLC = cPs/2 (see Figure 4 of Bogovalov 1999). Injection of energy at the termination intrabinary shock then proceeds with a periodicity of . The electrons accelerated to relativistic energies at different locations in the shock will radiate their energy on a timescale
where we used the relations for the synchrotron power emitted by a relativistic electron and to express the typical energy of synchrotron photons in terms of the electron energy . In the expressions above, is the Thomson scattering cross section and is the magnetic energy density. On the other hand, the light travel time between different regions of the shock is
As long as both the emission timescale and the light travel time between different locations of the emission regions are shorter than, say, half the spin period, a synchrotron-emitting spot (shaded in blue in Figure 14) will be seen to rotate coherently at the shock surface with a period PS/2. For a moderately large inclination angle, the exact value depending on the disk and shock relative height, the disk will absorb the emission coming from the spot closest to the observer (labeled as S2 in Figure 14). On the other hand, the emission from the spot located farther from the observer (S1) will be modulated sinusoidally as the spot rotates, as if the wind–disk shock were a sort of reflecting mirror. Two pulses of optical/X-ray synchrotron radiation will then be observed every spin cycle of the pulsar, the relative amplitude of which depends on the magnetic inclination angle, as well as on the viewing angle. Relativistic beaming and/or an ordered magnetic field would increase the anisotropy of the emitted radiation and the pulse amplitude. In this scenario, the large duty cycle of the observed X-ray pulse could result from the sum of the periodic emission emitted from different locations of the intrabinary shock, which are reached at different times by the spiraling-out current sheet, and are seen at different angles from the observer. Interestingly, the synchrotron timescale for optical photons is , compatible with the lag of optical pulses with respect to X-rays. The observed lag would find an immediate interpretation in our modeling, as optical synchrotron photons take a longer time to be emitted than X-rays.
The synchrotron timescale expressed by Equation (7) increases if the shock is located at a greater distance because the strength of the postshock magnetic field decreases linearly with distance (see Equation (6)). Eventually, it becomes comparable to half the spin period for 1 eV photons when the magnetic field at the shock is as low as 2 × 105 G. To produce optical coherent oscillations, the shock surface must then be located at k ≲ 2 times the light cylinder radius. The condition of the light travel time of different regions of the shock, , implies a similar constraint, . Remarkably, the latter condition is geometrical and does not depend on the energy of the photons. We speculate that the simultaneous disappearance of X-ray and optical pulses during the low mode might be due to the inner rings of the mass inflow being pushed outward by the pulsar wind, corresponding to an expansion of the termination radius beyond ≃π/(2 sin i) (k ≃ 2.2 for i = 45°) times the light cylinder radius. This is in agreement with the observation of radio flares during the X-ray low modes by Bogdanov et al. (2018), who interpreted them as episodes of ejection of optically thin plasmoids by the rotation-powered pulsar.
The synchrotron emission timescale expressed by Equation (7) is shorter than the flight time of accelerated particles in the shock region, rSH/c ≃ 33 (rSH/10 km) μs, provided that the latter is larger than a few kilometers. This ensures that the energy of the electrons is radiated away before they escape from the acceleration region.
The energy radiated in pulsed X-rays is 5 × 10−3 times the spin-down energy. The X-ray efficiency of typical isolated pulsar wind nebulae is usually of a few percent of the spin-down power (Kargaltsev & Pavlov 2010) and rarely reaches a value larger than 10% out of the reverberation phase (Younes et al. 2016; Torres & Lin 2018). Assuming an efficiency of the order of that observed from the Crab pulsar nebula (0.04), roughly 10% of the pulsar wind energy must be converted into X-rays to match the pulsed flux observed from PSR J1023+0038. For an isotropic distribution of the pulsar wind, the shock height must be hSH ∼ 0.2krLC to meet the energy requirement. This is 10 times larger than the height a Keplerian disk would have at a similar distance for a mass accretion rate , indicating that the shock must be vertically extended. An even more extended shock is required if the whole X-ray luminosity seen in the high mode () is due to synchrotron emission in the intrabinary shock, suggesting a higher efficiency of conversion of the pulsar wind in electron energy than previously assumed. A detailed modeling of the multiwavelength spectral energy distribution of PSR J1023+0038, including also the component observed at GeV energies, will be presented in a forthcoming paper.
In the context of a variable shock height hSH, an increase of the solid angle covered by the shock could explain the flares observed simultaneously in the optical and X-ray bands. As these sometimes reach a luminosity comparable to the spin-down power of PSR J1023+0038 (Bogdanov et al. 2015), it is evident that almost complete enshrouding of the pulsar by disk plasma and a conversion efficiency of spin-down power into electron energy close to unity would be required. It is unclear, however, why X-ray pulsations should disappear during flares. As we pointed out earlier, the nondetection could be due to the lower counting statistics in the X-ray band than in the optical band. The pulsed optical luminosity observed during flares is roughly a third of that observed in the high mode; a smaller decrease of the amplitude would be expected if flares were a sheer superimposition of unpulsed flux over the high-mode level. In addition, we observed a similar decrease in the optical pulse amplitude in intervals that formally fell in the high mode according to our definition, but were observed in between flares (i.e., toward the end of TNG2; see Figure 3). This would suggest that flaring intervals unlikely result from the addition of a component to the high-mode emission and should be treated separately. However, it is also possible that flares are produced in the outer regions of the disk, and the marked decrease of the optical pulse amplitude during flares could be due to the occurrence of low modes that cannot be identified from the X-ray light curve as they are outshined by the flaring emission. More observations of the pulsed amplitude decrease during flares are needed to break the degeneracy and identify the region where flares are produced.
The possibility that the accretion flow is stopped by the pulsar wind just beyond the light cylinder by the pulsar wind was recently explored with general relativistic MHD simulations by Parfrey & Tchekhovskoy (2017; see panel (d) in their Figure 4), who noted that X-ray emission would be expected from the trains of shocks and sound waves produced at the pulsar wind termination. They found that in this scenario, the amount of open magnetic flux was similar to the isolated pulsar case. The spin-down rate of the pulsar was then expected to be similar to that observed in the rotation-powered state. This is qualitatively in agreement with the ∼30% increase in the spin-down rate observed by J16 after the formation of an accretion disk, a factor much smaller than that expected if accretion and/or propeller ejection took place.
Ekşİ & Alpar (2005) showed that a stable equilibrium between the outward pressure of a rotation-powered pulsar and the inward pressure of the infalling matter can be realized if the termination shock is close to the light cylinder; this is due to the presence of a transition region from the near region inside the light cylinder, where the energy density of the electromagnetic field scales as , to the radiation zone far outside the light cylinder, where the ∼r−2 scaling holds. Equilibrium solutions for values of the disk truncation radius krLC with k ≳ 1 were found, with the exact value depending on the angle between the magnetic moment and the rotation axis, 1 − α. For instance, they found stable solutions with 1 < k < 2 for ξ = 10° and μ26 = 1, compatible with the assumptions of our model. On the other hand, far from the light cylinder, the radiation of a rotation-powered pulsar () decreases less steeply than the ram pressure of matter infalling under the gravitational pull of the compact object (). Because of this, a stable solution of a rotation-powered pulsar surrounded by an accretion disk would not ensue as the disk is expected to be fully ablated away by the pulsar wind as soon as rm ≫ rLC (Shvartsman 1970; Burderi et al. 2001). Takata et al. (2014) and Coti Zelati et al. (2014) modeled the multiwavelength emission of PSR J1023+0038 assuming that the rotation-powered pulsar wind interacted with the accretion disk to produce synchrotron X-ray emission far from the light cylinder (∼109 cm, i.e., k ∼ 125). Aside from the problems of stability that would probably ensue, we note that in the framework we propose coherent pulsations would not be produced at such a large distance from the pulsar.
Assuming that most of the pulsed emission is produced at the intrabinary shock between the pulsar wind and the disk does not rule out that magnetospheric rotation-powered pulses are still produced at a similar level than observed during the radio pulsar state. However, even if a radio pulsar were active, its pulses would be smeared and absorbed by material ejected by the pulsar wind (Stappers et al. 2014). On the other hand, the fivefold increase in gamma-ray flux observed after the disk formation could be due to Compton upscattering of disk UV photons off the cold relativistic pulsar wind, as proposed by Takata et al. (2014), while the pulsed magnetospheric emission keeps working at a similar rate to that in the radio pulsar state (Tam et al. 2010).
5. Conclusions
We presented the first simultaneous optical and X-ray high temporal resolution observations of PSR J1023+0038 in the accretion disk state, the only optical millisecond pulsar discovered so far (Ambrosino et al. 2017). We showed that optical pulsations are detected during the high mode observed in the X-ray light curve, in which X-ray pulsations also appear, while pulsations in both bands disappear in the low mode. Optical pulses are described by two harmonics, similar to X-ray pulses, and lag the X-ray pulses by ∼200 μs, although we caution that the absolute time calibration of SiFAP is based on just a single observations (see the Appendix). These findings suggest that the same phenomenon produces the pulses seen in both energy bands. Cyclotron emission from matter accreting onto the polar caps of the NS is not powerful enough to explain the pulsed optical luminosity. On the other hand, emission from the magnetosphere of a rotation-powered pulsar requires an unusually large efficiency of spin-down power conversion to match the optical and X-ray pulsed flux of PSR J1023+0038 with respect to other known pulsars. We argued that PSR J1023+0038 is a rotation-powered pulsar with a relativistic, highly magnetized wind interacting with the inflowing disk matter just beyond the light cylinder, creating a shock where the wind periodically deposits energy by accelerating electrons at the shock; these in turn produce optical and X-ray pulses through synchrotron emission. This would make PSR J1023+0038 the prototype of a few-hundred-kilometer-sized pulsar wind nebula and provide a unique opportunity to study the pulsar wind properties in the high magnetization regime rather than where they are particle dominated, as in the usual subparsec scale pulsar wind nebulae. This scenario also provides an explanation of the low mode and flares observed in the X-ray light curve in terms of the shock being pushed back by the pulsar radiation or increasing its size, respectively. Future observations will confirm the phase lag between optical and X-ray pulses, study the energy distribution of the pulses in the visible band, search for polarized pulsed emission, and look for gamma-ray pulsations, thus testing the scenario we proposed. On the other hand, magnetohydrodynamic simulations will be performed to demonstrate that pulsed emission can be indeed generated by a disk/wind intrabinary shock close to the light cylinder of a pulsar. The stability of a wind–disk shock just beyond the light cylinder over timescales of years should also be investigated. Other transitional millisecond pulsars in the subluminous disk state such as XSS J12270–4859 (de Martino et al. 2010, 2013; Bassa et al. 2014) and candidates like 1RXS J154439.4–112820 (Bogdanov & Halpern 2015), 1SXPS J042947.1–670320 (Strader et al. 2016), XMM J083850.4–282759 (Rea et al. 2017), and CXOU J110926.4–650224 (Coti Zelati et al. 2019), may be found in a similar state to that PSR J1023+0038. A search of optical pulsations in those source therefore seems warranted.
We are grateful to the XMM-Newton, TNG, and GTC directors for scheduling ToO observations in Director's Discretionary Time, as well as to all the teams involved for their effort in scheduling the simultaneous ToO observations of PSR J1023+0038. This work is partly based on observations made with the Italian Telescopio Nazionale Galileo (TNG), operated on the island of La Palma by the Fundación Galileo Galilei of the Istituto Nazionale di Astrofisica (INAF), and with the Gran Telescopio Canarias (GTC), which are installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, in the island of La Palma. Some of the scientific results reported in this study are based on observations obtained with XMM-Newton, which is a European Space Agency (ESA) science mission with instruments and contributions directly funded by ESA member states and NASA, with the Neil Gehrels Swift Observatory, which is a NASA/UK/ASI mission, with the NuSTAR mission, which is a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory and funded by NASA, with NICER, which is a NASA mission. Development of CIRCE was supported by the University of Florida and the National Science Foundation (grant AST-0352664), in collaboration with IUCAA. A.P. acknowledges funding from the EU's Horizon 2020 Framework Programme for Research and Innovation under the Marie Skodowska-Curie Individual Fellowship grant agreement 660657-TMSP-H2020-MSCA-IF-2014. We acknowledge financial support from the Italian Space Agency and National Institute for Astrophysics, ASI/INAF, under agreements ASI-INAF I/037/12/0 and ASI-INAF n.2017-14-H.0. L.S. acknowledges financial contribution from "iPeska" research grant (PI: Possenti) funded under the PRIN-INAF SKA/CTA presidential decree N. 70/2016. D.d.M. acknowledges support from "Towards SKA/CTA era" research grant (PI: M. Giroletti) funded under the PRIN-INAF SKA/CTA presidential decree N. 70/2016. A.P., F.C.Z., and D.T. acknowledge the International Space Science Institute (ISSI-Beijing), which funded and hosted the international team "Understanding and Unifying the Gamma-rays Emitting Scenarios in High Mass and Low Mass X-ray Binaries." A.P. acknowledges Cost Action PHAROS (CA16214), which funded the international workshop "The Challenge of Transitional Millisecond Pulsars," hosted by the INAF Osservatorio Astronomico di Roma. N.R. is supported by the ERC Consolidator Grant "MAGNESIA" (nr. 817661), and D.F.T., N.R., and F.C.Z., are supported by the Spanish Grant PGC2018-095512-B-I00, the Catalan grant SGR2017-1383, and COST Action "PHAROS" (CA 16124). A.V. was supported by grant 14.W03.31.0021 of the Ministry of Science and Higher Education of the Russian Federation. S.S.E. was supported in part by a University of Florida Research Foundation Professorship. Y.D. was supported in part by a University of Florida Graduate Student Fellowship. A.G. acknowledges L. Di Fabrizio and L. Riverol for the manufacturing of the original mask used in SiFAP. F.O. acknowledges the support of the H2020 Hemera program, grant agreement No. 730970.
Facilities: GTC (Circe) - , NICER - , NOT (ALFOSC) - , NuSTAR - , TNG(SiFAP) - , XMM. -
Appendix: SiFAP Observations of the Crab Pulsar
We observed the Crab pulsar with SiFAP at the 152 cm Cassini Telescope at Loiano Observatory for 3.3 ks starting on T0 = 57,724.04514 MJD. The target was observed at an average count rate of 45 × 103 s−1. We reduced the data following the same steps described in Section 2.6, correcting the arrival times for the assumed linear SiFAP clock drift, and reporting them to the solar system barycenter using the software Tempo2 (Hobbs et al. 2006), considering the position (R.A. 05h34m3197232 , decl. +22°00'52
0690 [J2000]) and the JPL DE200 ephemerides used by the Jodrell Bank monthly ephemerides (Lyne et al. 1993).21
We epoch-folded the optical time series in 1024 phase bins using the epoch of the nearest arrival time of the main pulse, TJB = 57,715.000000195995 MJD, the period 0.03372925321514(43) s, and the period derivative 4.197644(15) × 10−18. The maximum of the main pulse occurs at phase δϕCrab = −0.0054 ± 0.0005 (see Figure 15). Considering the uncertainty reported in the Jodrell Bank ephemerides to be 60 μs, we then estimate that the optical pulse lags the radio one by δτCrab = (181 ± 62) μs. This is compatible with the estimates given by Oosterbroek et al. (2008), who quoted (255 ± 21) μs from simultaneous optical and radio timing; Germanà et al. (2012; δτ ≃ 230 μs) Collins et al. (2012; δτ ∼ 178 μs); and Zampieri et al. (2014; δτ ∼ 240 μs). We then conclude that the absolute timing accuracy of SiFAP is better than 60 μs.
Figure 15. Pulse profile of the Crab pulsar observed with SiFAP at the Cassini Telescope at Loiano Observatory. Zero epoch corresponds to the nearest Jodrell Bank radio main pulse epoch TJB = 57,715.000000195995 MJD. The profile was obtained by folding the SiFAP time series in 1024 time bins using the period 0.03372925321514(43) s and the period derivative 4.197644(15) × 10−18 evaluated at TJB. The inset shows a magnification of the peak. The vertical dotted lines indicate the phase of the optical and radio maxima.
Download figure:
Standard image High-resolution imageA.1. The Optical Flux
To determine the optical flux of PSR J1023+0038 during the SiFAP/TNG observations, we first calibrated the flux/count rate conversion factor, k. The count rate observed by SiFAP is expressed as
where Fλ(λ) is the flux density, Ae = 9 × 104 cm2 is the effective area of the TNG mirror, R(λ) is the response of the detector equipped with a white filter (see Supplementary Figure 1 in Ambrosino et al. 2017), and E(λ) is the atmospheric extinction (we considered an air mass of 2.0 and 1.3 for TNG1 and TNG2, respectively). We measured the normalization factor C by evaluating Equation (9) for the reference star; we calculated the integral for the synthetic spectrum of a G2 V star (Pickles & Depagne 2010),22 normalized to give the observed flux density over the Sloan Digital Sky Service g filter (g = 15.86(1) mag, corresponding to 2.13 × 10−14 erg cm2 nm−1[=1.64 mJy]), and we matched it to the observed, background-subtracted count rate (RREF,TNG1 = 19456.8 s−1, RREF,TNG2 = 30750.2 s−1), obtaining CTNG1 = 0.434 and CTNG2 = 0.686. We attribute the difference between the normalization factors to the different atmospheric conditions over the two nights. The conversion factor k for PSR J1023+0038 was then estimated as
yielding k1 = 5.17 × 10−16 erg cm−2 and k2 = 3.25 × 10−16 erg cm−2 for TNG1 and TNG2, respectively. The rightmost column of Table 2 lists the 320–900 nm flux measured in the different modes by scaling the net observed count rates using these conversion factors and considering a 5% uncertainty. The dereddened flux in the 320–900 nm band is a factor of 1.22 larger than those values, considering the absorption column measured by Coti Zelati et al. (2014; NH = 5.2 × 1020 cm−2), a ratio AV/E(B − V) = 3.1, and the resulting color excess E(B − V) = 0.073 (Predehl & Schmitt 1995). We checked that performing this procedure to evaluate the ratio of the flux of PSR J1023+0038 to that of the reference star integrated over a B filter during the interval simultaneous to the TJO observation performed with the same filter (see Table 1), FJ1023/FREF = 0.43, is compatible with the flux ratio 0.41(2) indicated by the magnitudes observed (BPSR = 17.10(3) mag, BREF = 16.20(1) mag, FJ1023,TJO/FREF,TJO = 0.44 ± 0.02).
Footnotes
- 18
- 19
Available at http://www.eso.org/sci/software/eclipse/.
- 20
Note that Bednarek (2015) proposed a coexistence of an equatorial disk flow down to the NS surface and electron acceleration in high latitude slot gaps.
- 21
Available at http://www.jb.man.ac.uk/~pulsar/crab.html.
- 22
The spectrum is available at http://www.eso.org/sci/facilities/paranal/decommissioned/isaac/tools/.