Multiband Optical Light Curves of Black-widow Pulsars

"Black-widow" (BW) class binary pulsars lie at the extreme of the pulsar population. They have been recycled by mass transfer in a low-mass X-ray binary (LMXB) phase and are at present evaporating the companion remnant via a relativistic pulsar wind. For classical BWs, like the original PSR J1959+2048 system (Fruchter et al. 1988), the companion is a ∼0.02–0.05 M_⊙ substellar remnant of an LMXB mass donor in a P_b ≈ 2.5–10 hr orbit. In the closely related "Tidarren" subclass (Romani et al. 2016), the remnant is ≤0.02 M_⊙ and P_b < 2 hr, and the companion is the H-stripped core of a partly evolved star, evidently the descendents of the ultra-compact X-ray binaries (UCXBs). These systems are of great interest for the physics of millisecond pulsar (MSP) recycling and companion evaporation. Also, slow mass transfer during the X-ray binary phase allows substantial neutron star (NS) mass growth; some systems may even reach black hole collapse, so the heaviest BW avoiding this fate can probe the M_max NS limit. Because pulsar heating makes the substellar companion visible, optical observations can probe these phenomena. For example, optical spectroscopy can be used to measure the companion composition and radial velocity. Multicolor photometry, together with light-curve modeling, plays a key role in these studies since light-curve fitting constrains the distance, the companion Roche-lobe fill factor, and the heating. Most importantly, this fitting constrains the system inclination i, which together with spectroscopic measurements determines component masses ∝sin³i.

We have made measurements of a set of such MSPs using a variety of telescopes, including new observations, reanalysis of archival data sets, and use of published photometry. After summarizing the targets and the new data collected (Section 2), we describe our photometry (Section 3), which puts the measurements on a consistent (Sloan Digital Sky Survey (SDSS)-based) photometric system based on differential photometry against catalog stars in the field. We have focused on g-, r-, and i-band imaging, with some u- and z-band measurements. Additional images, especially optical, are on the Johnson system. Direct image photometry is supplemented in a number of cases with broadband magnitudes computed from differential spectroscopy. We generally have phase overlap between the direct and synthesized photometry, helping establish the zero-point (grayshift) of the synthesized points. After collecting the most complete multiband light curves possible, we fit (Section 4) the data with the ICARUS light-curve synthesis code (Breton et al. 2013), using the common assumption that the companions are directly heated by the pulsar radiation. In Sections 4.2–4.5 we discuss the results and their implications for companion heating. We conclude by considering some model amendments that might produce better agreement with theoretical expectations.

Searches in the directions of Fermi Large Area Telescope (LAT) γ-ray sources have discovered over forty "spider" (companion-evaporating) MSPs (BWs and their "redback" (RB) cousins with low-mass main-sequence secondaries). This is >10 times the pre-Fermi number, offering many opportunities for studying this dramatic phase in the life of interacting binaries. In Table 1 we summarize the basic observed properties of our nine systems. Figure 1 shows the BW/Tidarren population with our sample labeled in red.

Zoom In Zoom Out Reset image size

Table 1.  Black-widow Parameters

Name	Ps	Pb		x1	d	AV	μ	Range	fγ	Δγ
	(ms)	(hr)	(1034 erg s−1)	(lt-s)	(kpc)	(mag)	(mas yr−1)	i maga	b	b
J0023+0923	3.05	3.33	1.60	0.035	1.3	0.382	13.9	22.3-24.5	7.6	0
J0251+2606e	2.54	4.86	1.81	0.066	1.2	0.382	>14	22.8-	4.3	0.50
J0636+5128	2.87	1.60	0.58	0.009	1.1c	0.218	4.2	22.4-23.9	<1	0
J0952−0607	1.41	6.42	6.65d	0.063	1.7	0.137	⋯	22.0-24.4	2.2	0.20
J1124−3653	2.41	5.44	0.52	0.080	1.1	0.124	⋯	22.8-26.5	12.5	0.35
J1301+0833	1.84	6.48	6.65	0.078	1.2	0.027	26.9	21.4-24.8	7.6	0.40
J1959+2048	1.61	9.17	16.0	0.089	1.7	0.600	30.6	19.6-24.6	15.3	0.45
J2052+1219e	1.99	2.75	3.37	0.061	3.9	0.328	15.4	22.0-25.1	5.5	0.35
J2241−5236	2.19	3.50	2.50	0.026	1.0	0.040	17.2	21.3-24.0	26.6	0.55

Notes.

^aApproximate i magnitude range, when measured. For J1959, range in R. ^b4FGL Fermi LAT 0.1–100 GeV flux in units of double peak separation Δ_γ (0 for single peaks). ^cJ0636 timing parallax (Arzoumanian et al. 2018); all other distances are derived from the dispersion measure (DM; Yao et al. 2017). ^dFrom Nieder et al. (2019). ^eJ0251 and J2052 values from J. S. Deneva et al. (2019, in preparation). Other values are from atnf.csiro.au/research/pulsar/psrcat/ and astro.phys.wvu.edu/GalacticMSPs.txt.

Download table as:  ASCIITypeset image

Initial optical photometry of J0023+0923 (hereafter J0023) was published by Breton et al. (2013). PSR J0636+5128's (J0636) optical counterpart was described by Draghis & Romani (2018) and Kaplan et al. (2018). J0952−0607 (J0952) had some r imaging reported by Bassa et al. (2017). PSR J1301+0833's (J1301) counterpart was first described by Li et al. (2014) with additional photometry given by Romani et al. (2016). PSRs J0251+2606 (J0251) and J2052+1219 (J2052) were announced by Cromartie et al. (2016). J1124−3653 (J1124) was found in the GBT 350 GHz targeted survey (Hessels et al. 2011), while J2241−5236 (J2241) was discovered by Kieth et al. (2011). To the best of our knowledge the optical counterparts of these last four BWs have not been previously described, so in Figure 2 we present finder charts showing the companions near maximum brightness.

Zoom In Zoom Out Reset image size

Observations were obtained with a variety of telescope/camera systems. In Table 2 we summarize the data used in this paper. Starting from small aperture, on the Las Cumbres Observatory (LCO) 2 m telescopes we observed with the Sinistro cameras, obtaining gri measurements of our brighter sources near optical maximum. Unfortunately, these were near the system limiting magnitude; hence, low transparency, poor seeing, and tracking problems severely limited the utility of these data. In practice, LCO observations only contributed to the model fitting for J0023 and J2241 (Table 2), although other measurements were broadly consistent with more accurate photometry from other facilities. R-band photometry was obtained with the Tempelton camera on the MDM 2.4 m Hiltner telescope for J1301 on 2013 April 4 through 2014 May 28 (UT dates are used throughout this paper), as described by Li et al. (2014). J2052 was observed with the OSMOS spectograph in the imaging mode and the SDSS r filter on 2016 June 6.

Table 2.  Journal of Observations

Name	Tel./Inst.	Exp. (s)	MJD
J0023+0923	LCO2m	4 × 600 r	58366.057-.078
		11 × 600 r	58369.917-71.052
		1 × 600 g	58372.015
		3 × 600 g	58372.969-3.000
		4 × 600 g	58377.968-.997
	DEIMOS	2×600 R	57639.456-.463
	GMOSa	1 × 620 g, 2 × 320 i	55445.506-.518
		3 × 620 g, 4 × 320 i	55449.568-.617
	SOAR	11 × 600 i	58335.875-.946
		17 × 300 i	58336.869-.941
		2 × 300 g/r	58336.861-.876
J0251+2606	LRIS	200, 300 B/R,	57756.307-.939
		200, 6 × 300, 2 × 450 V,
		200, 3 × 300, 4 × 450 I
		13 × 320 g,	57991.566-.621
		6 × 300 R, 7 × 300 I
J0636+5128	GMOSa	10 × 420 g, 9 × 420 r	57012.324-.420
	NIRIa	21 × 60 K	56966.588-.610
		25 × 60 K	56999.422-.447
	NIRC2a	27 × 60 H, 28 × 60	56342.354-.408
	LRISb	32 × 450 g, r, i, z	58454,58455
J0952−0607	LRIS	19 × 310 g, 12 × 300 I	58455.501-.576
		5 × 300 R
	LRISb	12 × 900 g, r, i, z	58454.526-.639
	LRISb	4 × 900 g, r, i, z	58577.392-.426
J1124−3653	GMOS	3 × 300 r, 2 × 300 i	57417.197-.228
		2 × 300 r, 1 × 300 i	57420.260-.290
		8 × 600 g, 8 × 300 r, 7 × 300 i	57429.213-.323
		3 × 423 r	56727.365-.377
		11 × 300 i	57870.151-.241
J1301+0833	LRIS	2 × 300 I	57515.490-.509
	LRISb	21 × 600 g, r, i, z	57427-57572
	GMOS	3 × 420 r	56793.969-.985
	SDSS	1 × 54 g, r, i	55730.391
J1959+2048	R07	B, V, R, I, K	see R07
J2052+1219	LRIS	2 × 300 g/I	57305.264-.305
		3 × 300 g, 4 × 300 I	57640.346-.366
		2 × 300 u/Rs	57640.346-.366
		6 × 320 g,7 × 300 I	57991.46-.99
	DEIMOS	2 × 300 R, 3 × 300 I	57686.236-.249
	MDM	23 × 300 r + 2 × 900 r	57545.343-.469
	SAMI	30 × 180 i	58336.73-.80
J2241−5236	LCO 2 m	9 × 600 g, r, i	58310.127-4.260
	SOI	6 × 20, 14 × 180 i	57307.596-.603
	GHTS	3 × 180 g, r, i,	56514.728-5.851
		4 × 300 r, 1 × 300 i
	GMOS	3 × 600 r, 3 × 300 i	58365.740-.769

Notes.

^aRemeasured archival data. ^bSpectral synthesized magnitudes; see text.

Download table as:  ASCIITypeset image

Southern Astrophysical Research Telescope (SOAR) observations in 2017 were carried out with the Goodman High Throughput Spectrograph (GHTS) in imaging mode. In 2018, we observed with the SOAR Adaptive Module Imager (SAMI) which includes a UV laser system for correction of the ground-layer seeing (GLAO). A substantial portion of the run had poor seeing, where this GLAO system was of limited use, but several periods of subarcsecond seeing allow the corrections to produce a full width at half-maximum intensity (FWHM) as small as 05. Gemini GMOS data were obtained both from the science archive (for J0023, J0636, and J1301) and from new observations (under programs GS-2015B-FT3, GS-2017A-FT10, and GS-2018B-FT102).

Keck I/Low Resolution Imaging Spectrograph (LRIS) observations allowed us to get simultaneous red and blue images (largely g/I, but some V/I and u/R_s). A few Keck II/DEIMOS R exposures were obtained at minimum light during other programs. We also obtained Keck I/LRIS long-slit spectroscopy of a number of BWs, during which we adjusted the slit position angle (PA) to place a nearby comparison star on the slit. After extraction and calibration of these spectra using an automated pipeline (Perley 2019), we can integrate across broadband filters using the Image Reduction and Analysis Facility (IRAF) sbands script to obtain simultaneous relative magnitudes for the comparison star and the BW companion. These differential magnitudes can be calibrated with the (constant) catalog magnitudes of the comparison star, to correct for variable seeing and slit losses. (However, if the slit PA is not perfectly correct, there will be a grayshift, since the slit losses will be different for the two targets.) This gives simultaneous broadband magnitudes, even when the spectral signal-to-noise ratio (S/N) was quite low.

We obtained such spectral-synthesis magnitudes in the griz bands for J0636 (Δm = +0.2 mag), J1301 (Δm = 0.1 to 0.5 mag), and J0952 (Δm = 0.65, − 0.55 mag). The Δm gray offsets above were determined for each epoch and applied to the resulting BW magnitudes. We estimated uncertainties from the scatter in the resulting magnitudes near maximum brightness: J0636 (g, 0.06; r, 0.03; i, 0.05; z, 0.06), J0952 (g, 0.05; r, 0.02; i, 0.03; z, 0.06), and J1301 (g, 0.1; r, 0.04; i, 0.03; z, 0.06 at maximum brightness, doubling near quadrature). The epoch-dependent gray-offset estimates above are also uncertain by as much as Δm ≈ 0.1 mag. These spectral extractions are flagged in Table 2. At the brighter phases the exposures can have sufficient S/N for true spectral studies; these will be reported elsewhere.

In general, after standard charge-coupled device processing (bias subtraction and flat-fielding), we used IRAF to measure the FWHM of a number of bright stars to determine image quality. Then, identifying a grid of ∼15–50 well-detected stars in each frame, we registered the frames to subpixel accuracy. After measuring the position of the pulsar companion in several frames near optical maximum where it was best detected, we performed forced point-spread function (PSF) weighted photometry at the (fixed) position of the pulsar and our bright stars. The Gaussian function for the PSF weights was set to width 1.3 × FWHM for each frame and flux was extracted from a 1.7 × FWHM aperture. We then determined the catalog magnitudes of the grid stars in the relevant filter, using SDSS magnitudes directly when possible. When not, we used PanSTARRS g, r, i, z magnitudes, converting to SDSS using the transformations recommended by Finkbeiner et al. (2016):

where the a_i coefficient vectors are as follows: u (0.04438, −2.26095, −0.13387, 0.27099) [correcting from g_P1], g (−0.01808, −0.13595, 0.01941, −0.00183), r (−0.01836, −0.03577, 0.02612, −0.00558), i (0.01170, −0.00400, 0.00066, −0.00058), and z (0.08924, −0.20878, 0.10360, −0.02441). For two of our southern pulsars (J1124 and J2241) we have neither SDSS nor PanSTARRS coverage, and so we were forced to use AAVSO Photometric All-Sky Survey (APASS) photometry to transfer from bright stars to our infield star grid. We obtained g, r, i frames with the LCO 1 m system, where the large field provided more APASS stars for calibration of the nearby fainter comparison star grid.

These BWs have magnitude ranges as large as δm = 5 mag (Table 1). At maximum brightness, the photometric statistical errors can be unrealistically small, <0.01 mag. Therefore, we add a systematic σ_s = 0.03 mag in quadrature to all uncertainty estimates (except for those of the brightest BW, J1959) before the light-curve fitting to account for unmodeled zero-point, flat-fielding, and aperture errors. Conversely, at minimum brightness some objects, while formally detected in our forced (fixed position) astrometry, have very large statistical uncertainties (σ_m ≈ 1 mag). While we do include these points in our model fits, all detections with σ_m > 0.5 mag should probably be properly viewed as upper limits on the companion flux.

After barycentering the time of all observation midpoints and computing the orbital phases with the latest ephemerides available, we plot the folded light curves in Figures 3–5, with two cycles shown. Phase 0 is at the pulsar ascending node, so that pulsar inferior conjunction (near optical maximum) is at ϕ_B = 0.75.

We fit these data with the ICARUS modeling code (Breton et al. 2012) in its original form, where a bolometric luminosity L_X is assumed to propagate linearly from the pulsar point source and heat the companion. However, most of the pulsar spindown power is emitted as a relativistic magnetized e^±/B wind, rather than as photons. This wind is typically reprocessed in an intrabinary shock, and the resulting energetic particles can irradiate the surface or even be ducted by companion magnetic fields to surface caps. Such effects are described by Romani & Sanchez (2016) and Sanchez & Romani (2017), but for simplicity and consistency we here focus on the basic direct-heating fits; such model extensions may, however, be important for more realistic fits.

The direct-heating model takes as its basic parameters L_X, the source distance d, the companion Roche-lobe fill factor f_c, the uniform underlying temperature of the companion T_N, and the binary inclination i. L_X is reported assuming isotropic emission over the sky, but we suggest a model-dependent beaming correction for this in Section 4.1; the heating effect is from the flux directed at the companion. Our fit source distance can be compared with the radio dispersion measure (DM) distance estimate from Table 1. The fill factor is the ratio of the companion radius at the nose to the L₁ radius. Note that for the substellar BW companions there is no nuclear energy source. However, tidal effects can provide some internal heat and conduction must carry some energy from the day side to the night side, so T_N should be small but nonzero. In the model fitting, we use the precise orbital period P_b and the projected pulsar semimajor axis x₁ from the radio and/or gamma-ray ephemeris. For J1301 we can use the published estimate for the companion radial-velocity amplitude (Romani et al. 2016). Similarly, for J1959, van Kerkwijk et al. (2011) give , where the observed velocity is corrected by a geometrical correction factor K_cor ≈ 1.1 to account for the difference between the radial amplitude of the spectrally weighted center-of-light motion and the desired center-of-mass motion. With this we get an estimate of the mass ratio q = KP/(2πx₁).

For our other BWs, no radial velocity is yet published; we must assume an NS mass. Antoniadis et al. (2016) found that the MSP distribution is bimodal with a dominant peak at M_NS ≈ 1.4 M_⊙ and a subpopulation with M_NS ≈ 1.8 M_⊙ (width ). Strader et al. (2019) find that essentially all RB MSPs were from the higher mass peak. Here we assume a round M_NS = 1.5 M_⊙ to acknowledge that evolution and accretion probably increase the BW masses above the normal pulsar 1.4 M_⊙, without tying to a specific value. If one assumes that BWs, like RBs, populate the higher mass peak, our companion masses would increase by ∼13%. A more serious estimate of the pulsar mass distribution should be attempted once more BW radial velocities are available. This requires careful treatment of the radial velocity correction factor, so we do not discuss it further here.

The observed colors are affected by the intervening extinction A_V. For many of our objects, we can adopt the extinction estimates from the three-dimensional dust model of Green et al. (2015) available at argonaut.skymaps.info, using the Bayestar2019 model and A_V = 2.73 × E(g−r) values at the estimated distance. For the two southern sources, we can use as an upper limit the Schlafly & Finkbeiner (2011) total extinction column through the Galaxy from the NASA Extragalactic Database.

To obtain observed-band fluxes from the companion temperature distribution, we use a set of model-atmosphere spectra folded through the responses of our broadband filters, as provided by the Spanish Virtual Observatory (svo.cab.inta-csis.es). As in the past, we use the BTSettl atmospheres of Allard et al. (2012), to cover the very large temperature range from the heated day side to the cool night side, assembling a table of F_k(T, g) surface fluxes in each band k. Since spectroscopic studies of BW (and RB) companions indicate that the stripping can leave remnants with appreciable metal enrichment, we also explore the metallicity dependence of the model fits, using the available log [Z] = −1, −0.5, 0, +0.3, and +0.5 BTSettl models. Finally, we identified a flux-normalization error in our previous version of the code that resulted in d estimates being too large in Sanchez & Romani (2017); this has been amended and the direct-heating distances derived here should supersede those values.

As noted above, we add σ_s = 0.03 mag in quadrature to each point before fitting (except for J1959), to reduce sensitivity to the very small statistical photometric uncertainties at maximum brightness. We start the fitting with the Python SciPy least-squares minimization and then map the minimum and covariances with Markov Chain Monte Carlo runs (MultiNest v3.10). Characteristic corner plots from the latter are shown in Figure 6. The model filter light curves for the resulting best fits are shown in Figures 3–5 while the fit parameters are listed in Table 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  Binary Model Fits

Name	i	fc	LX/1033	TN	d	mc	χ2/DoF	fi/f*		ηheu	ηγ	ηH
	(deg)		(erg s−1)	(K)	(kpc)	(a)		(b)	(c)
J0023+0923							335/62(5.41)	0.46/0.48	1.15	0.29	0.18	0.12
J0251+2606							72/32(2.25)	0.74/0.81	0.86	0.34	0.47	0.23
J0636+5128							557/166(3.36)		0.56	0.42	0.02	0.10
J0952−0607							116/70(1.66)	1.16/0.89	6.65	0.12	0.07	0.27
J1124−3653							56/40(1.40)	1.16/0.99	0.52	0.44	2.48	0.54
J1301+0833							161/126(1.28)	1.00/0.89	2.53	0.20	0.15	0.18
J1959+2048							140/89(1.57)	0.78/0.79	8.93	0.11	0.07	0.27
J2052+1219							121/84(1.44)	0.73/0.79	1.11	0.30	0.68	0.68
J2241−5236							163/62(2.64)	0.87/0.86	1.77	0.24	0.24	0.14

Notes.

^aCompanion mass in 10⁻²M_⊙, as inferred from the fit assuming M_NS = 1.5 M_⊙ (except for J1301, J1959; see text). ^b"Outer gap" (OG) conversion between the Earth line of sight (f_i, at inclination i) and equatorial plane (f_*) isotropic flux estimates and the corresponding true sky average. When f_i/f_* > 1, the model predicts more companion heating than Earth-directed gamma-ray flux. ^cFor pulsars with proper motion μ this spindown power (units of 10³⁴ erg s⁻¹) is Shklovskii-corrected for I₄₅ = 1 and our fit distance; see Section 4.3.

Download table as:  ASCIITypeset image

The overall matches to the light curves are reasonable, with strongly heated companions brightest near ϕ_B ≈ 0.75, the subpulsar points. However, the χ²/Degrees of Freedom (DoFs) are relatively large for many of these fits. This may indicate that the photometric uncertainties are underestimated or that there is random variation about the quiescent light curve (e.g., flaring as seen for PSR J1311−3430; Romani et al. 2015). It is also likely that there are physical effects in the companion heating not captured in the basic model.

Five of our BWs (J0251, J0636, J1124, J1301, and J1959) have systematically decreased χ²/DoF using atmosphere tables with higher metallicities. Two (J0952 and J2052) appear to prefer log [Z] ≈ 0 or lower. Interestingly, these two BWs have hottest day-side temperatures in our sample, so it may be that in this higher T range the colors are less sensitive to metal overabundance. The last two objects show no significant dependence. The models only extend to log [Z] = +0.5, or ∼3× solar, but for the strongest dependence, in J1124, we can extrapolate from the most metal-rich models to infer a minimum at ∼4.5× solar abundance, although all one can really conclude is that log [Z] > 0.5 should produce even better fits for several of our pulsars (Figure 7).

Zoom In Zoom Out Reset image size

4.1. Spindown Power and Beaming Corrections

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 as (independent of binary properties). This means that the GeV efficiency is ∼0.01–0.3 for most γ-ray pulsars.

The spindown power should be corrected for the Shklovskii effect

where is the nominal spindown luminosity in units of for a moment of inertia , μ_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 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" for the sky output of γ-rays and 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.

4.2. Individual Data and Fit Results

4.3. Comparison of Efficiencies and Beaming Corrections

4.4. Inclination Distribution

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°.

4.5. Companion Masses and Radii

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

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.

We have collected a new set of photometric observations of BW MSPs and assembled multicolor light curves. We fit these data with a basic radiative heating model. On the whole, the fits allow such direct pulsar irradiation, but there are some challenges to the model. First, in a number of cases the reduced χ² is large. We suspect that underestimation of photometric errors is largely responsible, but a few objects (e.g., J0023) show evidence of nonthermal emission and others (e.g., J0251, J0952) may be variable. Such effects are not part of this basic model. The implied radiative heating efficiencies of 10%–60% are somewhat larger than assumed by other authors in the past.

We also compare the required photon heating with the observed LAT GeV flux, which dominates the pulsar photon emission along Earth's line of sight. The estimated γ-ray efficiencies η_γ are expected to be comparable to the "heuristic" γ-ray efficiency η_heu, and to the (beaming-corrected) heating efficiency η_H for a direct-heating model. The measured values do not always agree. For J0636 we can attribute the small η_γ to strong beaming away from Earth. But for J1124, the very large inferred η_γ probably indicates an overestimate of d. Conversely, for J0952 and J1959, the fit companion heating is more than twice the observed GeV luminosity at the photometric fit distance. While the inferred GeV flux would increase with d, this would rapidly exceed the expected heuristic gamma-ray luminosity. So for these objects a larger beaming correction (gamma-rays more tightly beamed at the companion) seems to be a more likely solution. Interestingly, the observed γ-ray pulse widths suggest larger i for individual objects in this model. This may indicate some bias in the direct-heating i fits, but it would be useful to also compare with the predictions of other gamma-ray beaming models.

The fits also give interesting information on the companions, which substantially fill their Roche lobes, with companion densities moderately larger than the minimum from their orbital period. Some objects were fit with substantially larger distance than the DM estimates (J0251, J1124). As noted above, these require unrealistically large γ-ray fluxes (unless the γ-rays are tightly beamed toward Earth, as can happen at larger i). Conversely, if the DM distances are adopted, the photometry demands small companion radii and large companion densities for these two objects, an evolutionary challenge.

The light-curve shapes also hint that other physical processes contribute to the optical emission. Small phase shifts in the optical maximum can be seen (e.g., J0636), and a few have maxima broader than the best-fit model curves in some bands (e.g., J1124) or have blue color minima deeper than the models (e.g., J2052). All of these suggest that some additional heat source or heat transport mechanism may be present. Since the distribution of measured i values is somewhat weighted to intermediate latitudes, and since the γ-ray pulse shapes generally prefer i larger than found in our fits, we may suspect that these effects introduce some bias into the measurements. It is important to continue to search for such effects as decreased i result in increased M_NS measurements for these BWs. As these are among the heaviest NSs known, this is important for equation-of-state constraints.

Overall, our photometry shows that the observed GeV flux is a plausible source of the direct heating of BW companions. But this is true only if the photometric distance estimates are substantially wrong in several cases. Since the fit distance also affects the companion size and density estimates and, less directly, the inferred system inclination, it is important to eliminate distance errors. We can of course continue to address the source of companion heating with improved multicolor photometry and spectroscopy, but it seems as if robust independent distances, ideally from pulsar timing or very long baseline interferometry parallax measurements, will be an important tool in understanding BW heating.

We thank D. Kandel, who helped with the ICARUS code, and the referee, whose careful reading and detailed comments helped us improve the paper. We are also grateful to Julia Deneva and Thankful Cromartie, who shared updated ephemerides for J0251 and J2052 in advance of publication. Daniel Perley and Brad Cenko kindly allowed us to use (and assisted with) their Keck LRIS imaging and spectroscopy reduction pipelines.

We are grateful to the staff at the Keck Observatory for their assistance. The W. M. Keck Observatory is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration (NASA); it was made possible by the generous financial support of the W. M. Keck Foundation. We extend special gratitude to those of Hawaiian ancestry on whose sacred mountain we are privileged to be guests. This work is based in part on observations obtained at the MDM Observatory, operated by Dartmouth College, Columbia University, Ohio State University, Ohio University, and the University of Michigan.

R.W.R. was supported in part by NASA grant 80NSSC17K0024. The research of A.V.F.s group was supported by the Christopher R. Redlich Fund, the TABASGO Foundation, and the Miller Institute for Basic Research in Science (U.C. Berkeley).

Software: ICARUS (Breton et al. 2012; Romani & Sanchez 2016), IRAF, SciPy (Jones et al 2001).