High-energy Variability of PSR J1311-3430

We first performed a gamma-ray spectral analysis. We used a $R=10^\circ$ region of interest (RoI), and analyzed the data following the binned likelihood analysis as described in Fermi Science Support Center (FSSC).⁶ When fitting the data, we used a power law with exponential cutoff model (PLEXP), ${dN}/{dE}={N}_{0}{(E/{E}_{0})}^{-{{\rm{\Gamma }}}_{1}}\exp (-{(E/{E}_{{\rm{c}}})}^{{{\rm{\Gamma }}}_{2}})$, with ${{\rm{\Gamma }}}_{2}$ held fixed at 1 for J1311 and fit all the bright sources (≳5σ) in the aperture, the Galactic diffuse emission (gll_iem_v06.fits; Acero et al. 2016), and the isotropic diffuse emission (iso_P8R2_SOURCE_V6_v06.txt; Ackermann et al. 2015) in the 100 MeV–300 GeV band with pylikelihood. In the fit, we also include all the 3FGL sources within 20° (Acero et al. 2015). We consider energy dispersion in the analysis,⁷ and therefore processed the data in the wider 60 MeV–500 GeV band. The best-fit parameters for J1311 are ${{\rm{\Gamma }}}_{1}=1.87\pm 0.03\pm 0.04$, ${E}_{{\rm{c}}}=5.3\pm 0.5\,\pm 0.5\,\mathrm{GeV}$, with 100 MeV–300 GeV flux of $9.2\pm 0.3\,\pm 0.2\times {10}^{-8}\ \mathrm{phs}\ {\mathrm{cm}}^{-2}\ {{\rm{s}}}^{-1}$ (see Table 1), consistent with the previous LAT spectral analysis within the uncertainties (Acero et al. 2015). The first and second uncertainties given for each parameter represent 1σ statistical and systematic uncertainties, respectively. Note that the systematic uncertainties are estimated by varying the scale for the response functions⁸ and normalization of the interstellar diffuse model (gll_iem_v06) by 6%. We then calculated changes in the best-fit parameters due to each systematic component and summed the errors in quadrature. We also varied the model; we hold all the parameters except for those of J1311 fixed at the 3FGL values or varied the number of sources to fit. We find that the results are all consistent with one another within the uncertainties. We further tried different RoI's, $R=5^\circ$ and $R=15^\circ$, and performed similar studies. The results for these new RoI's are also consistent with those for $R=10^\circ$ within the uncertainties.

Table 1.  Fit Results for the Orbital-phase-averaged and Orbital-phase-revolved Fermi-LAT Data in the "OFF" Pulse-phase Interval

Spin-phase-resolved and Orbital-phase-averaged
Phase	${{\rm{\Gamma }}}_{1}$	${E}_{{\rm{c}}}$	Fa
(GeV)
ALL	1.87 ± 0.03 ± 0.04	5.3 ± 0.5 ± 0.5	9.2 ± 0.3 ± 0.2
P1	1.71 ± 0.05 ± 0.06	4.7 ± 0.7 ± 0.4	14.5 ± 0.6 ± 1.4
Inter	1.67 ± 0.07 ± 0.11	3.7 ± 0.7 ± 0.6	8.1 ± 0.5 ± 1.3
P2	1.80 ± 0.04 ± 0.05	5.4 ± 0.8 ± 0.4	18.6 ± 0.6 ± 1.6
OFF	2.24 ± 0.12 ± 0.20	2.0 ± 0.7 ± 0.4	4.5 ± 0.3 ± 1.3
Orbital-phase-resolved in the "OFF" pulse-phase interval
Dip	2.54 ± 0.11 ± 0.11	3.3 ± 0.0b	4.3 ± 0.5 ± 1.2
Hump	2.33 ± 0.13 ± 0.19	3.3 ± 1.7 ± 0.6	5.2 ± 0.4 ± 1.3

Notes. 1σ statistical and systematic uncertainties are shown.

^a100 MeV–300 GeV flux in units of ${10}^{-8}\ \mathrm{phs}\,{\mathrm{cm}}^{-2}\ {{\rm{s}}}^{-1}$. ^bHeld fixed at the best-fit value for the Hump spectrum (see Figure 2).

Download table as:  ASCIITypeset image

To isolate the strong gamma-ray pulsations, we checked the timing solution of Pletsch et al. (2012), folding the Fermi-LAT events. These events were selected from small apertures (R ≲ 2°). The source probabilities of each event, based on energy and position, were computed using gtsrcprob, and the time series was phased and folded using tempo2. Although the pulsar spin modulation shows very strongly, this solution, based on earlier LAT data, leaves a large drift in the pulse arrival time residual especially after MJD 56600. We therefore derived a new solution covering all data until MJD 57510.

Starting from the timing solution of Pletsch et al. (2012), we used the PINT software package⁹ and an unbinned Markov chain Monte Carlo maximum likelihood sampling technique (Sanpa-Arsa 2016). We modeled the pulse profile using two asymmetric Lorentzians for the main peaks and a simple Gaussian for the bridge emission. The parameters of the template were varied jointly with the timing model parameters, allowing us to marginalize over the profile model. We found that accurately tracking the pulse phase required the addition of an orbital period derivative (${\dot{P}}_{B}$), proper motion (${\mu }_{\alpha }\cos \delta$ and ${\mu }_{\delta }$), and a second spin frequency derivative ($\ddot{\nu }$). Orbital eccentricity is undetected and fixed at e = 0. The resulting best-fit parameters and their 1σ confidence intervals from the posterior distribution are reported in Table 2, and the resulting pulse profile appears in Figure 1.

Zoom In Zoom Out Reset image size

Table 2.  Timing Parameters for PSR J1311−3430

R.A. (J2000)	${13}^{{\rm{h}}}{11}^{{\rm{m}}}45\buildrel{\rm{s}}\over{.} 72358(1)$
Decl. (J2000)	$-34^\circ 30^{\prime} 30\buildrel{\prime\prime}\over{.} 342(3)$
${\mu }_{\alpha }\cos \delta$ (mas yr−1)	−6.8(6)
${\mu }_{\delta }$ (mas yr−1)	−3.5(8)
Position Epoch (MJD)	56228
ν (s−1)	390.56839299885(1)
$\dot{\nu }$ (s−2)	$-3.1882(3)\times {10}^{-15}$
$\ddot{\nu }$ (s−3)	$1.03(8)\times {10}^{-25}$
Epoch (MJD)	56228
Binary model	ELL1a
${P}_{{\rm{B}}}$ (day)	0.0651157347(2)
${\dot{P}}_{{\rm{B}}}$ (s s−1)	3.7(3) × 10−12
Projected semimajor axis (lt-s)	0.010575(2)
${T}_{\mathrm{ASC}}$ (MJD)	56009.129455(4)

Note.

^aSee Edwards et al. (2006) for the timing model definition.

Download table as:  ASCIITypeset image

By cutting on spin phase (${\phi }_{P}$) to isolate the off-pulse region, one can minimize the pulsar magnetospheric flux and search for nonmagnetospheric emission. In a study of 2500 days of Pass-7 LAT data Xing & Wang (2015) found evidence for an orbital modulation of the off-pulse $\gt 0.2\,\mathrm{GeV}$ flux, finding ∼3σ (0.003 chance probability) significance. To study such emission with our improved data set we constructed orbital light curves in the off-pulse (${\phi }_{P}=0.127-0.587$) and three on-pulse (P1, P2, and the interpulse, Figure 2) intervals. A concern for this study is that the orbital period of ${P}_{b}\simeq 94$ minutes is very close to that of the Fermi satellite (96 minutes), hence artificial modulation may be produced due to the satellite's orbital motion. We therefore folded time-resolved exposure (${\rm{\Delta }}T=30$ s) on the orbital period and confirmed that exposure variations across the phase bins used for the orbital light curves in Figure 2 were negligible (≲1%). We also confirmed that varying zenith angle cuts to adjust the imperfect exclusion of Earth-limb gamma-rays introduced no variation beyond the simple scaling due to changes in source count statistics.

Zoom In Zoom Out Reset image size

As expected, no on-pulse interval shows significant modulation. However, the off-pulse light curves show modulation consistent with that reported by Xing & Wang (2015), deserving further study. Because we use a source-weighted photon analysis, we did not need to use the very small apertures of that study. This avoids the extreme sensitivity of modulation significance to the aperture choice. With a conservative 2° aperture, we compute the cumulative weighted H-test significance (Kerr 2011) starting from the beginning or end of the data set; a significance that grows as more exposure is included indicates a real, steady signal. We track the modulation in low (200 MeV–1 GeV) and high (1 GeV–100 GeV) energy bands in addition to the full LAT range. The cumulative significance curves are shown in Figure 3, and the final low- and high-band light curves are shown in Figure 2. We further changed the aperture size (e.g., $R=3^\circ \mbox{--}5^\circ$) and found that the results do not change significantly. Although larger counts for larger apertures may result in higher significance, especially at lower energies with the probability weighting, the weighting may not be accurate at the event-by-event level for timing analysis, hence an improvement in significance is not necessarily expected. As noted above, the exposure variation rapidly becomes negligible, dropping to the percent level after ∼100 days of exposure. Monte Carlo simulations confirm that including this small variation in the bin exposure makes negligible difference to the H-test significance.

Zoom In Zoom Out Reset image size

We note that the nearby flaring source 3FGL J1316.0−3338 (hereafter J1316) at a distance of $d\sim 1^\circ$ may affect the results. Although it is not very likely that this source exhibits periodic modulation on J1311 periods, a short-time flare of J1316 may be a concern. This effect is minimized because we use probability weighting, but there may be some residual effects. We therefore produced a light curve for J1316¹⁰ and found two relatively large flares. We excised the time intervals for the these flares (±7 days), and recomputed the H-test significance for J1311 in the three energy bands as we did above. We find that the results do not change whether or not we excise the J1316 flare intervals.

The low-band signal fades at late times while the high-band signal continues the secular increase in significance; this suggests some unknown, slowly varying, soft contamination. The full-band signal shows the largest significance, reaching $p\lesssim 2\times {10}^{-4}$. This significance decreases slightly, with $\sim 4.5\times {10}^{-4}$ by the end of the data set; this is as expected since the low-band counts dominate. This $\sim 10\times$ improvement in significance over the Xing & Wang (2015) result is encouraging, but requires further exposure to be definitive. We can, of course, improve exclusion of low-energy and pulse contamination, by varying phase, aperture, time, and energy cuts. As can be seen in Figure 1, using a narrower phase interval may suppress contaminations of the pulse peaks better and hence reveal stronger modulation. We therefore varied the selection criteria including the phase interval, and find that the results are not very sensitive to the selection criteria. For example, using ${\phi }_{p}=0.182\mbox{--}0.549$, aperture $R=1\buildrel{\circ}\over{.} 27$, E_γ = 1.01–100 GeV and interval MJD = 56047–57431, we reach p ∼ 10⁻⁶ significance. Interestingly, this "tuned" high-energy light curve is nearly 100% modulated. However, given that the post-trials significance is $\gtrsim {10}^{-4}$, comparable to that for the initial untuned cuts, we cannot claim that such a large modulation is required.

For phase-resolved spectral analysis, we generated the phase-resolved spectra for four different pulse-phase intervals: off-pulse (0.127–0.587), P1 (0.597–0.737), interpulse (0.757–0.917), and P2 (0.937–1.107) intervals, and performed likelihood analysis using pylikelihood with the data of each pulse interval. For these, normalizations for all the sources in the 3FGL model file were scaled down by the phase interval and the resulting flux for J1311 was scaled up by the same amount. The results are rather insensitive to the precise definition of these phase intervals. Here, we held all the parameters (except for those of PLEXP for J1311) fixed at the best-fit values obtained in the phase-averaged analysis above. The results are shown in Table 1. We further tried to see if the result changes if we fit the diffuse backgrounds (gll_iem_v06.fits and iso_P8R2_SOURCE_V6_v06.txt) as well, and found that the results do not change significantly.

For the off-pulse interval, we fit both the total flux and the spectra from the "hump" phase (${\phi }_{B}=0.6\mbox{--}1.2$) and the lower "dip" phase (${\phi }_{B}=0.2\mbox{--}0.6$). The spectral parameters are not significantly different although the hump power law is nominally slightly harder (Table 1). Note that the cutoff energy ${E}_{{\rm{c}}}$ for the "dip" phase is not well constrained, so we hold ${E}_{{\rm{c}}}$ at the value obtained for the "hump" phase. The spectra are shown in Figure 4 and the fit parameters are presented in Table 1. We have also fit a simple power-law (PL) model to the off-pulse phase spectra. Using the log-likelihood fit statistic, we find that PLEXP fits better than PL does with ${\rm{\Delta }}\mathrm{log}{ \mathcal L }$ of 9, 1, and 7 for "all," "dip," and "hump," respectively, in the 100 MeV–300 GeV band. For the "all" and the "hump" spectra, a curvature (i.e., a cutoff) is still required (≳3σ), but we cannot tell whether the "dip" spectrum requires a curvature. Nevertheless, the shapes are very different from the on-pulse spectra (Figure 4).

Zoom In Zoom Out Reset image size

Because the low-energy light curves' significance appears concentrated in the early data, we also investigate long-term source variability. The variability index (Acero et al. 2015) in the 3FGL catalog is 52, and hence the total source flux has no significant variability. Variability might still be significant in some phase bins. As in Section 3.1, we selected four phase intervals and constructed a light curve for each phase interval, with time bin 2 Ms for the faint off-pulse interval and 1 Ms for the rest. Note that there is a nearby variable source (J1316) that can contaminate the J1311 light curves. For each time bin, we performed a likelihood analysis using pylikelihood to fit the amplitudes of J1311 and J1316. All other source and background parameters are held fixed at those obtained by a pulse-phase-resolved analysis (Section 3.1). We used the best-fit source fluxes to calculate ${\chi }^{2}$ for a constant flux. As expected the on-pulse phases dominated by magnetospheric emission are consistent with steady flux. The off-pulse phase gives ${\chi }^{2}/\mathrm{dof}=57/86$, so even this interval is consistent with steady emission. We conclude that there is no long-term time variability in J1311, which is consistent with the study of Torres et al. (2017).

J1311 is known to exhibit optical and X-ray flares (Kataoka et al. 2012; Romani 2012), and a possible orbital modulation in the X-ray band has been reported (Kataoka et al. 2012; Arumugasamy et al. 2015). So a re-examination of the X-ray observations provides an important comparison with the gamma-ray modulation.

We re-examined archival X-ray observations taken with Swift, XMM-Newton, and Suzaku. Source fluxes were extracted from apertures with radii $R=40^{\prime\prime}$, $R=16^{\prime\prime}$, and $R=40^{\prime\prime}$, respectively, while the background was monitored by nearby source-free apertures of radii $R=60^{\prime\prime}$, $R=32^{\prime\prime}$, and $R=60^{\prime\prime}$. Event times were barycentered with the ephemeris of Section 3.1. In Figure 5, we plot the binned (${\rm{\Delta }}t=300$ s) XMM-Newton and Suzaku flux time series, with t = 0 set at the ascending node prior to observation start. We use Chandra observations for spectral analysis only because of the short exposures and small number of counts. The source and background events are extracted using a $R=3^{\prime\prime}$ aperture and an annular region with ${R}_{\mathrm{in}}=5^{\prime\prime}$ and ${R}_{\mathrm{out}}=10^{\prime\prime}$, respectively.

Zoom In Zoom Out Reset image size

We investigate source variability in the time series using the Bayesian-block algorithm (Scargle 1998; Scargle et al. 2013) implemented in the python astroML package (Vanderplas et al. 2012). The algorithm computes the number of optimal blocks and the change points for the time series. In the source time series, we find 15, 13, 5, 7, and 2 blocks for XMM-Newton, Suzaku/AO6, Suzaku/AO4, Swift, and Chandra; there are some blocks with much larger count rates than the minimum level and others similar to the minimum level (Figure 5). We performed the same analysis with the background time series to see if some of the source variabilities can be attributed to background activities. The Suzaku/AO4, Swift, and Chandra backgrounds are well explained with a single block (i.e., no variability), and the others require two to four blocks, suggesting some variability in the background; these variabilities are small and do not seem to correlate with the source activity. Furthermore, the background flux is only a small fraction of the source flux, hence the small background variabilities are unlikely to drive the source activities.

Strong flaring variability is clearly visible (Figure 5), with episodes reaching $\gt 10\times$ the quiescent flux. We have marked some flux levels to help guide the eye at 4× the quiescent flux (red horizontal lines). There is no obvious preferred phase for flaring events. It appears, especially in the XMM-Newton data, that the flares can be clustered in episodes lasting several orbits of ∼10–20 ks.

The Swift data are snapshots over many years, so all we can show is the phase dependence of intervals (100–1400 s) where the flux was $\gt 4\times$ the quiescence (Figure 5(e)). In Figure 5(e), for only a few intervals was this excess more than 90% significant (red points in Figure 5(e)). For the significance, we calculate the Poisson probabilities of having the observed counts in the time bins greater than the background plus $4\times$ the (XMM-estimated) quiescent counts, considering the trial factor (number of data bins in the figure). We denote these time bins with green lines in the Bayesian-block light curve as well (Figure 5(f)). Two of the high-significance data points in Figure 5(e) have a corresponding peak in Figure 5(f). The lowest one does not have a counterpart (green line at T ≈ 7 × 10⁴ s in Figure 5(f)); this block has slightly elevated flux compared to the minimum, suggesting low level variability. In Figure 5(f), the first flare does not have a high-significance (red) counterpart in Figure 5(e); the corresponding point is denoted as a green cross in Figure 5(e). Again there is no clear phase grouping of the limited number of significant flares.

XMM-Newton has the highest count rate and the longest continued coverage. To test for orbital modulation, we formed the light curve for "flare" ($T\lt 56872.1$ MJD, the first ∼4 orbits in Figure 6 left) and "quiescent" ($T\gt 56872.1$ MJD) periods (Figure 6 left). The quiescence light curve is fully consistent with constant flux (p = 0.5); no variability is seen in our Bayesian-block analysis. With limited flare events in XMM alone the light curve is spiky, but we cannot discern a preferred phase. Substantially longer integrations are needed to infer the detailed flare behavior—but the data already show that flares can start at any phase.

Zoom In Zoom Out Reset image size

The XMM-Newton data provide by far the best statistics for measuring the X-ray spectrum. We can distinguish between the "flare" epoch and "quiescent" epochs. For each, we extracted source and background events using apertures of $R=16^{\prime\prime}$ and $R=32^{\prime\prime}$, respectively. We then calculated the response files using the standard tools of SAS version 20141104_1833 (rmfgen and arfgen). We then jointly fit the two 0.3–10 keV spectra with an absorbed power-law model having a common absorbing column density (${N}_{{\rm{H}}}$). The source absorption is undetectably small; it is consistent with 0 with the current statistics. Unsurprisingly, the flare spectrum is substantially harder than that in quiescence (Table 3).

Table 3.  Absorbed Power-law Fits to the X-Ray Data

State	${N}_{{\rm{H}}}$	Γ	${F}_{0.3\mbox{--}10\mathrm{keV}}$
	(${10}^{20}{\mathrm{cm}}^{-2}$)	⋯	(${10}^{-14}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$)
XMM flare	1 ± 2a	1.26 ± 0.07	47 ± 2
XMM quiescence	1 ± 2a	1.67 ± 0.09	6.6 ± 0.3
Suzaku AO4	1b	1.3 ± 0.1	26 ± 2
Suzaku AO4 Flare	1b	1.19 ± 0.15	94 ± 10
Suzaku AO4 Quie.	1b	1.6 ± 0.2	14 ± 2
Suzaku AO6	1b	1.16 ± 0.05	40 ± 2
Chandra 11790	1b	1.3 ± 0.3	18 ± 4
Chandra 13285	1b	1.4 ± 0.3	8 ± 2
Swift combined	1b	1.20 ± 0.15	19 ± 2

Notes.

^aTied. ^bFrozen.

Download table as:  ASCIITypeset image

While other X-ray data sets provide insufficient statistics for accurate fits, we can at least estimate fluxes and hardness. In all cases, we fix ${N}_{{\rm{H}}}=1\times {10}^{20}\,{\mathrm{cm}}^{-2}$. For the Suzaku data, we use $R=40^{\prime\prime}$ and $R=60^{\prime\prime}$ for the source and background apertures, respectively, and response files generated with XSELECT. The AO4 data had a clear bright flare and the separate flare and quiescence spectra have indices similar to the XMM measured values. In AO6, the source is relatively hard and bright when compared to the quiescent state observed with XMM-Newton, suggesting substantial flare contributions (Table 3; see also Figure 5(c)). In the first of two archival CXO ACIS data sets, the source is relatively bright with substantial flare contribution; the second is closer to quiescence. The spectral index uncertainties are too large to select a state. For Swift, we simply combined all XRT observations, used pre-processed files for the source spectra, and separately created background spectra using $R=60^{\prime\prime}$ apertures in the source-free regions. A fit to the merged spectrum gave an intermediate level average flux and a poorly determined, but hard Γ, again suggesting substantial contribution from poorly measured flares. Note that if the XMM flare/quiescence flux ratio is typical, a flare duty cycle as small as 10% will cause significant contamination in the spectral fits.

The optical light curve has the dramatic ∼4 mag ($\sim 40\times$) orbital modulation characteristic of BW pulsars, with maximum at pulsar inferior conjunction ${\phi }_{B}=0.75$ when the heated face of the companion is best visible. In addition, it shows optical flares with rise times as short as ∼300 s and amplitudes as large as $\sim 6\times$ the peak flux (Romani 2012). We would like to know how these optical events relate to the X-ray flares above.

We are fortunate to have simultaneous optical data from the Swift/UVOT and XMM-Newton/Optical monitor (OM). Unfortunately, these data were taken with a variety of integration times and filters. The XMM-Newton/OM exposures were 2 ks (${\rm{\Delta }}{\phi }_{B}=0.36$). This coarse time sampling makes it difficult to identify any but the largest flares. The Swift/UVOT exposure varied between 30 s and 1.6 ks. To extract fluxes for an initial study of the optical behavior, we used standard tools (uvotsource and omichain; Section 2).  For the Swift data, we barycenter correct the exposure midpoint times. For the XMM-Newton data, we used the barycenter-corrected X-ray photon arrival times, and interpolate them to obtain the barycenter-corrected OM times. The omichain automatically detects the source in the XMM-Newton/OM data and properly calculates the source flux. For the Swift/UVOT data, we used a $R=5^{\prime\prime}$ and a $R=15^{\prime\prime}$ apertures for the source and the background, respectively, to measure the source flux.

The resulting barycentered midpoints of the Swift observations are phased with the binary ephemeris and plotted in Figure 6, with detections, mostly near ${\phi }_{B}=0.75$, in the first cycle and upper limits in the second cycle. To know if these detections are bright compared with the quiescent state, we need to compare with the quiescent (pulsar heating) light curve. We can use the ground-based observations of Romani et al. (2015a) to interpolate quiescent V and B light curves. However, the quiescent UV fluxes of J1311 are not known. The hot ${T}_{\mathrm{eff}}\approx {\rm{12,000}}$ K face of the He companion must contribute some flux. However, after dereddening by an estimated $E(B-V)=0.173$ (Galactic average toward J1311; Romani et al. 2015a), we find that the phase average optical UV spectrum does not match a simple blackbody. This is in contradiction to Kataoka et al.'s (2012) UVOT analysis.

Still, we need an estimate of the quiescent flux in the UV filters, and so, to be conservative, we scaled the best-fit model B light curve assuming a simple blackbody spectrum across the UVOT bands to generate U, UVW1, UVM2, and UVW2 light curves; these estimates are shown in Figure 6 and are also consistent with the Swift upper limits. Note that the modest Swift V sensitivity means that the source was only detected in a high (flare) state. One B detection also likely represents a mild flare. In the UVM2 and UVW2 bands, the measured fluxes are, in general, higher than the blackbody-extrapolated model fluxes. This may suggest that the source spectrum at high frequencies deviates from a simple blackbody. Several UV detections during the "night" phase (${\phi }_{B}=0.0\mbox{--}0.5$) clearly represent a large increase over the expected quiescent flux. The XMM-Newton/OM data are very coarse and have few detections, yet we can identify several bright flare points here, as well.

We have identified flare intervals in each X-ray data set (${\rm{\Delta }}t\,\sim$ 300–400 s) for which the flux is $\gtrsim 4\times$ the quiescent value; adjacent flare bins are counted as a single flare. When finding the number of flares, we use the light curves made with constant time bins instead of the Bayesian-block representation because we are concerned with the gap removal; multiple flares separated by data gaps may be seen as one when the gap is removed. Note again that "flare" means the time intervals with measured flux greater than $4\times$ the quiescent level in our definition. We note that this definition for flare is arbitrary; we use this 4× the quiescent level only for reference flux (e.g., for "large" flux), and hence the "flare" definition should be considered with care. We then histogram the start phase of these flare events as a function of orbital phase (Figure 7). While it is clear that flares can be visible at any phase, with this coarse binning, we see a mild, ∼50% enhancement at and just before ${\phi }_{B}=0.75$, pulsar inferior conjunction when we view the heated face most directly. However, this excess has low (p = 0.2) significance, and so additional monitoring will be needed to determine the flare phase distribution. Note that the difference in exposure time for each phase bin in Figure 7 is ∼10% (minimum to maximum), and has no significant impact. We show the exposure corrected flare counts in Figure 7.

Zoom In Zoom Out Reset image size

We would like to know the duty cycle of this flaring activity. For this, we use the Bayesian-block representation of the light curves (Figures 5(b)–(d) and (f)). In the light curves, we calculate the count rates and the exposure times for each time bin. If we adopt the $\gt 4\times$ flare definition, then fully 19% of all intervals are in a flaring state. However, many of these must be statistically insignificant fluctuations. If we additionally require a 90% significance for the interval to have a flux excess using the Poisson statistic with the background and the XMM-estimated quiescent source flux, the flare duty cycle drops to 7%. Note that the duty cycle inferred this way is similar to what we infer using the binned light curves. We conclude that the true duty cycle at this $4\times$ threshold lies between these generous (19%) and conservative (7%) rates.

We also show the number of flares as a function of flux. In Figure 7 (bottom), we plot the distribution of flare fluxes from all observations, with the flare flux in units of the quiescent flux for the given observatory. Of course some excesses are low significance, and to account for this we also plot the flare flux weighted by the inverse of the flux uncertainties. This better approximates the true distribution, which is evidently a rather steep ($N\sim {\mathrm{fl}}^{-3}$) power law.

Furthermore, we can investigate the association of X-ray and optical flares of J1311. Because XMM-Newton and Swift provide simultaneous optical observations, data taken with these telescopes are particularly useful for this study. The OM data of XMM-Newton were taken only with the V filter during the flare period, and we find that J1311 was detected only when the X-ray flaring activity is strongest. The OM time bins are of course large, but during the early bright flare V-band detections seem to span all phases, including "night" phases. We thus infer that these V detections are flare associated. The other OM filters were used during quiescence. The detections in these bands are only during the day phase, suggesting that we see only the heating modulation of the companion.

The many Swift observations allow us to plot the excess above background X-ray flux against the excess optical/UV flux beyond the heating model (Figure 8). The figure shows a correlation between X-ray and optical fluxes. The significance of the linear correlation is 7σ (although this is, of course, sensitive to the UV flux model that defines the flare excess). If we separately consider each band, the significance is highest in the UVM2 and the UVW2 bands, being $4\sigma$ and 5σ, respectively. We believe this suffices to show that the X-ray and optical/UV flares are due to the same events. Recently, Halpern et al. (2017) noticed correlated X-ray and optical variability in BW candidate XMMU J083850.38−282756.8, which may represent similar flaring events.

Zoom In Zoom Out Reset image size

Finally, we would like to know if these X-ray/optical flare events are associated with the off-pulse GeV modulation shown in Figure 2. The low-significance modulation in Figure 7 peaks at the rising part of the off-pulse LAT light curve, so such association is at least plausible. However, we are unable to associate individual flare events with gamma-ray off-pulse activity. For $R\lt 1^\circ$ from the source position, energy >1 GeV, and pulse phase intervals of 0.127–0.587, we get a total count rate of 0.17 day⁻¹, including background. However, the integrated flare time from all X-ray observations to date is 1.1 day. None of these gamma-rays arrive during an X-ray flare, but the Poisson probability of this is 80%. In fact, given the substantial duty cycle for the flares, even if all off-pulse LAT photons arose during flare periods, we would expect ∼1 coincident count. We attempted to boost the off-pulse count rate, by opening the energy range to ≳100 MeV (1.4 off-pulse photons/day), but this gives only 1.6 expected coincident counts (with a substantial fraction being background events). Two events are actually seen, no strong excess. Our only inference is that the Fermi-LAT flux during the X-ray flare intervals is ≲10× the quiescent flux (90% upper limit). Much more extensive (30–100 day) X-ray or optical flare monitoring is needed to test for an event-level gamma-ray association. Thus we cannot make a direct connection and, for the foreseeable future, only an indirect connection of the associated phase seems possible.