Multiband Optical Light Curves of Black-widow Pulsars

Because we assume that the companion is directly (radiatively) heated by the pulsar, it is useful to compare the results of the photometric heating model with other measures of the pulsar output. In fact, the νF_ν spectral energy distributions of rotation-powered pulsars are strongly dominated by the GeV peak observed by the Fermi LAT, so the observed LAT fluxes (Table 1) provide an estimate of the incident photon flux (scaled with the distance). Heuristically, the LAT GeV luminosity scales with pulsar spindown power $\dot{E}$ as ${L}_{\mathrm{heu}}\approx {({10}^{33}\mathrm{erg}{{\rm{s}}}^{-1}\dot{E})}^{1/2}$ (independent of binary properties). This means that the GeV efficiency ${\eta }_{\mathrm{heu}}\approx {({10}^{33}\mathrm{erg}{{\rm{s}}}^{-1}/\dot{E})}^{1/2}$ is ∼0.01–0.3 for most γ-ray pulsars.

The spindown power should be corrected for the Shklovskii effect

$$\begin{eqnarray*}&&{\dot{E}}_{c}=({10}^{34}\,\mathrm{erg}\,{{\rm{s}}}^{-1}){\dot{E}}_{34}{I}_{45}[1-0.0096{\mu }_{\mathrm{mas}\,{\mathrm{yr}}^{-1}}^{2}{d}_{\mathrm{kpc}}/({\dot{E}}_{34}{P}_{\mathrm{ms}}^{2})],\end{eqnarray*}$$

where ${\dot{E}}_{34}$ is the nominal spindown luminosity in units of ${10}^{34}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ for a moment of inertia ${10}^{45}{I}_{45}\,{\rm{g}}\,{\mathrm{cm}}^{2}$, μ_mas/yr is the total proper motion in mas yr⁻¹, P_ms is the spin period in ms (from Table 1), and d_kpc is the fit (or DM) distance. This correction can be quite substantial; for example, at the nominal DM distance J2052's spindown power would be reduced by a factor of 3, decreasing the energy budget available for companion heating.

As the model fits estimate this heating luminosity, it is useful to compare with the heuristic expectation L_heu and the observed gamma-ray power. However, the LAT gamma-ray flux is that directed at Earth's line of sight at inclination i, while the heating flux is directed at the companion star, which for spin-aligned MSPs will be at the pulsar equatorial plane. Since in most models pulsar emission is not isotropic, we must correct these fluxes to isotropic all-sky equivalents with "beaming corrections" f_i and f_*, respectively, to make such comparisons.

We expect the pulsar magnetosphere flux to be directed toward the spin equator, as seen in numerical models (e.g., Tchekhovskoy et al. 2016). There the energy/momentum flux of the pulsar wind is concentrated as $\propto {\sin }^{n}{\theta }_{* }$ with n = 2 or even n = 4 for an oblique rotator (θ_* is the angle from the pulsar spin axis, assumed perpendicular to the orbital plane). This flux is largely particles and fields, but the direct (photon) radiation should be similarly concentrated. Gamma-ray pulsar models typically have the emission occur near the light cylinder or even in the near-field wave zone outside the light cylinder (e.g., Kalapotharakos et al. 2018). In such models, the radiation beamed from these active zones is also thought to be equatorially concentrated. Note that there may also be (evidently fainter) emission beamed along the dipole magnetic axis, because lower-altitude radio pulses imply pair production (and hence beamed gamma-ray emission) in this zone, as well.

For a concrete example, we use the geometrical "outer gap" (OG) model whose beaming factors f_Ω(α, ζ, η_heu) = f/f_iso are given by Figure 16 of Romani & Watters (2010) as a function of magnetic inclination α, viewing angle ζ = i, and heuristic γ-ray efficiency η_heu. These are meant to correct observed fluxes to full-sky isotropic emission equivalences, but they also provide the model-dependant correction for the flux viewed at i to the flux directed at the companion near ζ = π/2 (see Table 3, right section). Our fit i has an uncertainty range (and the companion subtends nonzero angle), so the correction f_i is computed by weighting over the i uncertainty range and scanning over the α range which (for that i) allows Earth to see the γ-ray beam. For the companion heating we compute f_* by integrating over the angle subtended by the companion, and scanning the allowed α. We assume that angles are distributed isotropically within the allowed ranges. Two examples are shown in Figure 8 and the values are listed in Table 3. Note that these values are averaged over a range of angles and so, understandably, are rather close to 1; true values for a specific pulsar's α and i can be quite far from unity.

Zoom In Zoom Out Reset image size

We find it convenient to divide the (beaming-corrected) gamma-ray flux and heating flux by the (Shklovskii-corrected) spindown power and thus refer to the gamma-ray "efficiency" ${\eta }_{\gamma }={f}_{\gamma }{f}_{i}4\pi {d}^{2}/{\dot{E}}_{c}$ for the sky output of γ-rays and ${\eta }_{H}\,={L}_{x}{f}_{* }/{\dot{E}}_{c}$ for the output as estimated from the heating flux (see Table 3). If the measurements were perfect and the OG model for the GeV beaming were correct, these would be the same. As η_γ is an Earth-based estimate it depends on d, while η_H is measured in the BW system and so is (after Shklovskii correction) distance-independent. Inconsistency indicates bias in the light-curve fits, imperfections in the beaming model, or both. We evaluate the consistency in Section 4.3.

PSR J0023+0923 had a few GMOS images reported by Breton et al. (2013), which we remeasured and used in our fits. We made a number of LCO 2 m queue observations, but tracking in this region was exceptionally poor and so only a small fraction of the exposures provided usable photometry. A few Keck/DEIMOS exposures provided higher S/N detections in R, and two long sequences of SAMI i imaging provided good coverage at maximum brightness. These i data are interesting as an abrupt step occurs near maximum at ϕ_B = 0.8; this feature is consistently present on two consecutive nights in the middle of continuous image sequences. It is hard to imagine such sharp light-curve features from any thermal surface emission; one possible origin is beamed emission from an intrabinary shock. Additional observations, especially simultaneous multicolor photometry and sensitive X-ray light-curve measurements, can test whether there is a contribution from nonthermal emission at this phase.

The few points near minimum brightness (especially the GMOS g values) drive our solution to large inclination i. The new i-band magnitudes support this solution, but imply a relatively large night-side temperature. Breton et al. (2013) noted that this BW does not display radio eclipses and with a handful of points they find a smaller i = 58° ± 14°, which is ∼1.5σ away from our fit minimum, requiring Δχ² = +60 with our data. The J0023 fit has the worst reduced χ² of our set, apparently due to an under-prediction of the g maximum, and an indication in g and r that the maximum is narrower in the bluer colors. Additional photometry at minimum, perhaps in the z band, would help constrain T_N. Higher S/N multicolor light curves are needed to map the maximum and constrain the nonthermal emission, if present.

PSR J0251+2606 has poor light-curve coverage and, lacking detections near minimum brightness, the inclination is not well constrained. Moreover, the relatively good photometry near maximum suggests significant short-term variability, especially in the bluer V and G data. The R and I light curves are relatively well behaved. The most puzzling aspect of the fit is the large distance (2.8 times the DM estimate), and the resulting large inferred γ-ray luminosity and direct heating power, which are 47% and 23% of the spindown power, respectively (see Section 4.3). At present we do not have a full proper motion, but Shklovskii correction using the lower limit on the right-ascension component already substantially increases the implied efficiency. The γ and heating flux estimates would be in agreement at d ≈ 2.3 kpc, so a smaller distance seems likely. Clearly additional measurements, especially at minimum brightness, are needed.

PSR J0636+5128 has a relatively shallow light curve, and both Draghis & Romani (2018) and Kaplan et al. (2018) inferred a small inclination i, as also found here. An early timing parallax (Stovall et al. 2014) gave a very small distance, but this has been superseded by more extended timing, giving $d={1.1}_{-0.3}^{+0.6}$ kpc (Arzoumanian et al. 2018), and our photometric distance fit is in good agreement with this value. While not in the 4FGL source catalog (hence the upper limit in Table 1), it does have a LAT pulse detection (Smith et al. 2019) and is thus a (faint) GeV emitter. The required direct heating L_X is only 10% of the spindown power, but exceeds the observed γ-ray flux, as discussed below. Note that with this very small i, the estimated companion mass is not especially small, with ${M}_{c}\approx 0.018{({M}_{p}/1.5{M}_{\odot })}^{2/3}$. If the pulsar is massive, J0636 is closer to the BW range than most Tidarrens. Spectroscopic measurement of the companion composition and radial velocity are needed to test whether it is the Tiddaren descendent of an UCXB. The large χ²/DoF is mostly due to a significant phase shift (data late by Δϕ = 0.02; applying this arbitrary shift reduces χ² by ∼200).

PSR J0952−0607: In Bassa et al. (2017), a partial r Isaac Newton Telescope (INT) light curve showed the expected ϕ_B = 0.75 maximum, but placed few constraints on system parameters. We obtained (2018 December) Keck/LRIS spectra across maximum brightness that provide griz colors, supplemented by two-color LRIS imaging during minimum (on the next night), and a few more LRIS spectra in 2019 April. These provide appreciable color data. We determined the spectral gray shifts by a match to the INT r light curve, showing that the slit placement was not always optimal. The last three imaging points (2 × gI, 1 × gR, marked in Figure 3) are substantially brighter than any reasonable extrapolation of the prior day's light curve—we infer that J0952 was undergoing a flare event at ϕ_B ≈ 0.4–0.5, such as seen in other BWs. Also, relatively poor sky subtraction beyond 8500 Å, affected the z-band estimates of the last two spectral measurements, so we exclude these points. As this paper was being submitted for publication, Nieder et al. (2019) presented observations of J0952 that include ULTRACAM/HiPERCAM photometry. Although presented in flux rather than magnitudes, it appears that these data also include a similar (red) flaring event, in a similar phase range.

J0952's peak appears relatively narrow; thus, our fits prefer small f_c with minimal peak-broadening due to ellipsoidal modulation. This is the only object for which the fits prefer f_c < 0.5; for such small values the companion is nearly spherical and the distance is highly covariant with f_c, with little constraint on either. The very best fit occurs for log[Z] = −0.5, at which the solution drives to unphysically small d and f_c. This appears to be due to a z-band offset of ∼ −0.02 mag for this metallicity, which better matches the data. For both smaller and larger log[Z] the fits give d ≈ 3–4 kpc, but the d minima are shallow, allowing smaller distances down to the DM estimate of 1.7 kpc with Δχ² ≈ +8. However, larger d are relatively strongly excluded; with the 5.6 kpc of Nieder et al. (2019) we find Δχ² = +50. We conclude that the very small d of the log[Z] = −0.5 fit may be an artifact of imperfect photometric calibration and adopt here the agnostic (but slightly worse) log[Z] = 0 fits. Future precision photometry, especially in the near-infrared (or a parallax measurement), will be needed to study this important object. All parameters other than d and f_c were nearly independent of metallicity, so we consider these to be relatively robust.

PSR J1124−3653 was observed in several visits through the Gemini Fast Turnaround Queue program. One observational challenge is the very bright R = 12 mag star at the same decl., whose saturated core caused a charge-transfer-efficiency trail to pass very close to the pulsar in initial observations (Figure 2). For exposures capturing minimum brightness, we placed this star behind the guider pick-off mirror to avoid this issue. Our ICARUS fit is quite good, although the g data peak appears somewhat broader than that of the model. The main challenge is a fit distance ∼2.5 times larger than the DM estimate. As the spindown luminosity is small (even before Shklovskii correction), the required γ-ray heating power is 2.4× the nominal spindown luminosity. The γ-ray flux and direct-heating power at d = 1.2 kpc are in better agreement with the DM estimate. However, a direct-heating fit with d fixed at this value is much worse, with χ²/DoF ≈ 2.5 times larger, and i decreased in an (unsuccessful) attempt to fit the broad maximum. We suspect that a departure from simple direct heating is needed to accommodate a smaller distance.

PSR J1301+0833's light curve has points converted to r from Li et al. (2014), but the rest are the imaging and synthesized magnitudes discussed by Romani et al. (2016). Here our fit is not very sensitive to log [Z], but we do achieve a lower χ² and a closer distance estimate than found in the latter paper (although still larger than the DM estimate). At this 2 kpc distance the corrected spindown power is ${\dot{E}}_{c}=2.5\,\times {10}^{34}\,{I}_{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, so the heating flux requires an acceptable 18%. The fill factor remains reasonable and so this solution should supersede that of Romani et al. (2016), although the implied 260 km s⁻¹ space velocity is still unusually large for an MSP.

PSR J1959+2048 is the original BW MSP. Its pre-Fermi discovery likely stems from its relatively large P_b = 9.17 hr, because substantial time is thus spent out of eclipse with an unobscured pulsar detectable at radio wavelengths. Interest in inclination fitting is especially high because the van Kerkwijk et al. (2011) measurement of a large companion radial velocity coupled with the Reynolds et al. (2007) estimate of i = 65° ± 2° indicates a very high M_p = 2.40 ± 0.12 M_⊙. We rely here on the optical/infrared photometry previously published by Reynolds et al. (2007). Refitting, we find a weak preference for high log [Z] and a revised distance of 2.3 kpc. The i is consistent with, but even slightly smaller than, the R07 estimate. At this distance the direct-heating flux is 27% of the corrected spindown power. Thus, like J1301, the revised fit is now consistent with direct heating, although the inferred heating power exceeds the γ-ray luminosity by a factor of ∼4. This strengthens van Kerkwijk et al. (2011)'s conclusion that, for direct heating, the NS mass is large, although this critical direct heating assumption may still be questioned, as we discuss in Sanchez & Romani (2017) and in future work.

PSR J2052+1219: The MDM statistical errors seem to underrepresent the true fluctuations, so these were doubled, dropping χ² from 214 to 121. Otherwise the fit is quite good. For consistency with the other objects, we did not fit the two u points, but they are close to the predicted model. The Keck r and i points are somewhat lower than the model, suggesting that the fit should have a higher i and/or a lower T_N. Interestingly, this BW has the strongest preference for low log [Z]. In this case if A_V is free it stays close to 0.3, independent of the model log [Z]. Spectra of this object may test the heavy-element abundance. The model fit implies a large distance quite consistent with the DM estimate. At this distance the inferred heating and γ-ray powers agree, but are a large fraction of the spindown luminosity.

PSR J2241−5236 is particularly interesting since An et al. (2018) found evidence for GeV orbital modulation, suggesting that an intrabinary shock contributes in this band. The GeV peak is single, implying a tangent view of the shock, as would occur for i ≈ 30°–50°, depending on the companion wind momentum. Our modest fit i is consistent with this interpretation.

Figure 8 shows that for most magnetic inclinations α we expect little or no GeV flux at viewing angles i < 30°. Thus, J0636, whose shallow light curve strongly supports the fit low inclination, provides a basic test of the beaming picture. Indeed, it is an anomalously weak γ-ray pulsar (Table 1), only detected since the LAT events could be folded on the pre-existing radio MSP ephemeris. The heuristic luminosity law predicts a γ-ray efficiency of η_heu = 0.42; the observed 4FGL γ-ray flux limit corresponds to η_obs < 0.02, so it is underluminous by a factor of ∼20. Since we do not expect any OG flux directed at our fit i = 23°, there is no meaningful f_i in Table 3. (The flux we do see may be fainter GeV emission from lower altitude, directed along the magnetic pole.) Note that the required heating efficiency is only 10% of η_heu, which is easily accommodated.

However, when we compare η_γ and η_H for our full sample, we find only a mild correlation (Spearman correlation cor = 0.35). Clearly, the beaming correction or the heating fit parameters are imperfect. Several factors can contribute toward this disagreement. The sizes of the efficiencies are an important clue. For example, for J1124 we have η_γ = 2.5. The moment of inertia may be larger than the standard I₄₅ = 1, which is plausible for heavy, highly recycled pulsars. But as noted above, an overestimate of the photometric distance likely dominates. Returning J1124 to its DM distance would decrease η_γ to a large but plausible 0.4, which is in good agreement with η_H. This example highlights the fact that comparison of η_γ and η_H, assuming an accurate beaming correction, gives an independent estimate of the distance. Such an estimate would decrease J0252's distance to ∼2.3 kpc, increase J0952's to ∼8 kpc, and increase J1959's to ∼4 kpc. However, it seems unlikely that distance errors are the full explanation. Beaming corrections, even if imperfect, are clearly needed for some BWs. Certainly J0636 requires such correction. Also, inspection of Figure 8 shows that for our fit inclination i = 45° for J0952, many α would not produce an OG GeV detection. Also, distance alone is unlikely to explain the disagreement for J1959, because at 4 kpc the implied γ-ray luminosity would be ∼2.5 times its expected heuristic level. Improved estimates of MSP beaming are therefore desirable.

In our present beaming estimates, we have been forced to average over α (horizontal green bands in Figure 8), which result in f_i and f_* close to unity, whereas corrections for particular angles can be quite substantial. Other observables might resolve this degeneracy. For example, the beaming model also predicts the phase separation Δ_γ of the principle γ-ray peaks. All of our BWs except J0023 and J0636 have a fairly well-defined double structure, so we could use this to constrain the α range (in our given inclination i range) most consistent with the γ-ray pulse shape. However, the large Δ_γ of most of our BWs is only obtained in this OG model for i > 60°. Thus, if this OG model is applicable, the i must be underestimated. Since this beaming model was computed for young pulsars rather than BWs, we might suspect that a pulse shape match will require a different model. Nevertheless, when applied to more sources and more modern MSP beaming maps, this direct view/radiative heating comparison provides a novel way to test γ-ray magnetosphere models.

We have seen that there is some tension between our fit i values and the expectations of the OG model. Figure 9 compares the photometric fit i values to a simple isotropic distribution and to the i expected when one observes within the OG beam of Figure 8. In the OG case the lack of very low i does boost somewhat the fraction expected to be viewed edge-on. The data seem to have an excess at i = 40°–60°, but these differences are not statistically significant, with the KS test allowing both models. More objects will be needed to chose between these hypotheses.

Zoom In Zoom Out Reset image size

If the mid-range i excess above proves significant, it may be that BWs are most easily discovered at such inclinations. For example, the bright GeV MSP emission may be beamed away from the spin (orbital) axis, but companion evaporation produces a dense ionized plasma wind in the orbital plane, resulting in eclipses, dispersion, and scattering that will make it more difficult to detect BW radio counterparts at i ≈ 90°. The competition between these effects could bias detections to i ≈ 40°–60°.

Taking our fit i at face value, we can replace the commonly used minimum companion masses with true masses. The resulting m_c, assuming a pulsar mass of 1.5 M_⊙, are plotted with their uncertainty ranges in Figures 1 and 10. The masses cluster around 0.03 M_⊙ and it seems that all BW companions are heavier than 0.017 M_⊙. The biggest shift is, of course, for J0636, whose revised M_c ≈ 0.018 M_⊙ encroaches on the classical BW range.

With our inclination and fill factor estimates we can examine companion densities for clues to the evolution. First we recall the minimum density for a Roche-lobe-filling companion

$$\begin{eqnarray*}&&\rho =0.129{P}_{b}^{-2}\,{\rm{g}}\,{\mathrm{cm}}^{-3},\end{eqnarray*}$$

with the orbital period in days. In Figure 10 we show these lower density limits for several orbital periods. For comparison we plot our photometric estimates of the masses and radii. Here the error bars for the mass uncertainties incorporate only the i measurement error and we assume M_NS = 1.5 M_⊙. The dashed lines connect to the M_NS = 2 M_⊙ position to display the mass-scale range. For J1301 and J1959 we can use the published companion radial velocity to set the mass scale. Since we wish to focus here on companion properties, we do not discuss the pulsar masses, but defer that to a future publication describing the spectroscopic modeling.

Zoom In Zoom Out Reset image size

Note that most companions are relatively close to Roche-lobe filling and hence near the minimum density for their orbital period. Displayed in this way, we see that the BWs in our sample seem to group into families; this appears as well in the P_B distribution, and with a larger sample we can probe the evolution that leads to these objects. Note also that some objects are uncomfortably close to their minimum density for an assumed M_NS = 1.5 M_⊙; this supports the conclusion of Strader et al. (2019) that RBs (and by extension BWs) are, as a class, massive NSs. J0952 is well within its Roche lobe at our best-fit photometric distance, and a direct-heating model has ρ ≈ 10 g cm⁻³, although with our photometric fill-factor/distance uncertainty the density may be as low as 6 g cm⁻³ or as high as 19 g cm⁻³. As noted above, the large best-fit distance for J1124 makes the inferred GeV efficiency unreasonably large at η_γ = 2.5. This can be mitigated by smaller distances, at a cost in decreased fill factor and increased companion density. At its 1.1 kpc DM distance we would find a companion density ∼97 g cm⁻³ (25× increase). This large density suggests evolution from Roche-lobe contact at P_B ≈ 0.9 hr, which seems implausible. Thus, while closer distances would help reduce the large η_γ, the distances are still likely to be larger than the DM estimate. A more substantial beaming correction f_i could also help bring the efficiency estimates into agreement.