Constraining Relativistic Bow Shock Properties in Rotation-Powered Millisecond Pulsar Binaries

TL;DR

This paper models the shock geometry and X-ray emission in millisecond pulsar binaries, constraining shock properties through radio and X-ray observations, and highlights the importance of Doppler boosting and shock stand-off distance.

Contribution

It introduces a geometric model linking radio eclipses and X-ray light curves to shock properties, providing new constraints on shock geometry in pulsar binaries.

Findings

Shock stand-off distance constrained by radio eclipses.

X-ray light curves consistent with Doppler boosting effects.

Model degeneracies suggest need for further transport physics.

Abstract

Multiwavelength followup of unidentified Fermi sources has vastly expanded the number of known galactic-field "black widow" and "redback" millisecond pulsar binaries. Focusing on their rotation-powered state, we interpret the radio to X-ray phenomenology in a consistent framework. We advocate the existence of two distinct modes differing in their intrabinary shock orientation, distinguished by the phase-centering of the double-peaked X-ray orbital modulation originating from mildly-relativistic Doppler boosting. By constructing a geometric model for radio eclipses, we constrain the shock geometry as functions of binary inclination and shock stand-off $R_0$. We develop synthetic X-ray synchrotron orbital light curves and explore the model parameter space allowed by radio eclipse constraints applied on archetypal systems B1957+20 and J1023+0038. For B1957+20, from radio eclipses the stand-off is $R_0 \sim 0.15$ -- $0.3$ fraction of binary separation from the companion center, depending on the orbit inclination. Constructed X-ray light curves for B1957+20 using these values are qualitatively consistent with those observed, and we find occultation of the shock by the companion as a minor influence, demanding significant Doppler factors to yield double peaks. For J1023+0038, radio eclipses imply $R_0 \lesssim 0.4$ while X-ray light curves suggest $0.1\lesssim R_0 \lesssim 0.3$ (from the pulsar). Degeneracies in the model parameter space encourage further development to include transport considerations. Generically, the spatial variation along the shock of the underlying electron power-law index should yield energy-dependence in the shape of light curves motivating future X-ray phase-resolved spectroscopic studies to probe the unknown physics of pulsar winds and relativistic shock acceleration therein.

Figures

Tables

Table 1. Table 1: Rotation-Powered BWs and RBs with DP X-Ray Light Curve Morphology

Name	Type	DP Phase Centering	Refs.
B1957+20	BW	SC	[1]
J0024–7204W	RB	IC	[2]
J1023+0038	RB	IC	[3]
J12270–4859	RB	IC	[4]
J1723–2837	RB	IC	[5]
J2039–5618	RB	IC	[6]
J2129–0429	RB	IC	[7]
J2215+5135	RB	IC	[8]
J2339.6–0532	BW	IC	[9]

Equations

$$B_{\mathrm{w}}\approx {\left(\frac{3{\dot{E}}_{\mathrm{SD}}}{2c}\right)}^{1/2}\frac{1}{r_{\mathrm{s}}}=22{\left(\frac{{\dot{E}}_{\mathrm{SD}}}{10^{35}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{1/2}\left(\frac{10^{11}\mathrm{cm}}{r_{\mathrm{s}}}\right)G.$$

$$⟨{\gamma }_{\mathrm{w}}⟩\gtrsim {\left(\frac{3{\dot{E}}_{\mathrm{SD}}}{2c}\right)}^{1/2}\frac{e}{{\mathcal{M}}_{\pm }2m_e{c}^2}\approx \frac{7\times 10^8}{{\mathcal{M}}_{\pm }}{\left(\frac{{\dot{E}}_{\mathrm{SD}}}{10^{35}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{1/2},$$

$$B_{\mathrm{s}}\gtrsim B_{\mathrm{s},\min }\approx 4.4{ϵ}_{X,\max }^{1/3}{\left(\frac{10^9\mathrm{cm}}{r_{\mathrm{L}}}\right)}^{2/3}\mathrm{G},$$

$$R_{\mathrm{iso}}(\theta )$$

$${\theta }_{\infty }-\tan {\theta }_{\infty }=\frac{\pi }{1-{\eta }_{\mathrm{w}}},$$

$$R_0=\frac{\sqrt{{\eta }_{\mathrm{w}}}}{1+\sqrt{{\eta }_{\mathrm{w}}}}.$$

$${\eta }_{\mathrm{w}}$$

$$\xi =\frac{1}{4\pi }{\left(\frac{R_0}{1-R_0}\right)}^2\delta \Omega -\frac{2\pi R_{\ast }^2{\sigma }_B{T}_{\mathrm{cold}}^4}{{\dot{E}}_{\mathrm{SD}}}\approx \frac{1}{4\pi }{\left(\frac{R_0}{1-R_0}\right)}^2\delta \Omega ,$$

$$B_{\ast }\lesssim 4\times 10^3{\left(\frac{{\dot{E}}_{\mathrm{SD}}}{10^{35}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{1/2}\left(\frac{10^{11}\mathrm{cm}}{a}\right){\left(\frac{0.05}{R_{\ast }/a}\right)}^3{\left(\frac{R_0}{0.25}\right)}^3\left(\frac{0.75}{1-R_0}\right)\mathrm{G},$$

$$\cos {\theta }_{\max ,\mathrm{R}}=\cos i\cos {\vartheta }_{\mathrm{sb}}+\sin i\sin {\vartheta }_{\mathrm{sb}}\cos (\pi f_E).$$

$$f_E=\frac{1}{\pi }\arccos \left(\frac{\cos {\theta }_{\max ,\mathrm{R}}}{\sin i}\right)$$

$$p_{\Gamma }\equiv \Gamma {\beta }_{\Gamma }={\left(\Gamma \beta \right)}_{\max }\left(\frac{\theta }{{\theta }_{\max ,\mathrm{X}}}\right),$$

$$\beta (\theta )=\frac{{\left(\Gamma \beta \right)}_{\max }\theta }{\sqrt{{\left(\Gamma \beta \right)}_{\max }^2{\theta }^2+{\theta }_{\max ,\mathrm{X}}^2}},\Gamma (\theta )=\frac{1}{\sqrt{1-\beta (\theta {)}^2}}.$$

$${\tau }_{\mathcal{D},\perp }\gtrsim \frac{{\ell }^2}{2{\kappa }_{\parallel }}=1.4\times 10^4{\left(\frac{\ell }{10^9\mathrm{cm}}\right)}^2\left(\frac{B_{\mathrm{s}}}{1\mathrm{G}}\right){\left(\frac{v_{\mathrm{p}}}{10^9\mathrm{cm}{\mathrm{s}}^{-1}}\right)}^{-2}\mathrm{s},$$

$${\sigma }_{\mathrm{shocked}}\lesssim \frac{B_{\mathrm{s},\min }^2}{4\pi n_{\mathrm{e},\mathrm{MSP}}⟨{\gamma }_e⟩m_e{c}^2}\approx \frac{3\times 10^6}{{\mathcal{M}}_{\pm }⟨{\gamma }_e⟩}{\left(\frac{B_{\mathrm{s},\min }}{1\mathrm{G}}\right)}^2{\left(\frac{r_{\mathrm{s}}}{10^{11}\mathrm{cm}}\right)}^2{\left(\frac{10^{35}\mathrm{erg}{\mathrm{s}}^{-1}}{{\dot{E}}_{\mathrm{SD}}}\right)}^{1/2}.$$

$${\bar{j}}_{{ϵ}_f}=\frac{d\bar{E}}{d\bar{V}d\bar{t}d{\overline{ϵ}}_{f}d{\bar{\Omega }}_f}=m_e{c}^2{\overline{ϵ}}_f{\bar{\dot{n}}}_{\mathrm{SR}}^{\mathrm{iso}},$$

$${\bar{j}}_{{\overline{ϵ}}_f}\approx \mathcal{C}(p){\left(\frac{3}{2}\right)}^{(p+1)/2}\frac{{\Sigma }_{\mathrm{e}}(\theta ,\varphi )}{\Delta (\theta ,\varphi )}{\sigma }_{T}c{U}_{\mathrm{cr}}{\left(\frac{{\bar{B}}_s(\theta ,\varphi )}{B_{\mathrm{cr}}}\right)}^{(p+1)/2}{\overline{ϵ}}_f^{-(p-1)/2},$$

$$\mathcal{C}(p)=\frac{2^{(p-1)/2}\sqrt{3}}{32\sqrt{\pi }}{\left[\Gamma \left(\frac{p}{4}+\frac{7}{4}\right)\right]}^{-1}\Gamma \left(\frac{p}{4}-\frac{1}{12}\right)\Gamma \left(\frac{p}{4}+\frac{19}{12}\right)\Gamma \left(\frac{p}{4}+\frac{1}{4}\right).$$

$$L_{\mathrm{SR}}({\Omega }_f,{ϵ}_f)\equiv \frac{d{\mathcal{L}}_{\mathrm{SR}}}{d{ϵ}_{f}d{\Omega }_f}=\int dV{ϵ}_f^2\left(\frac{j_{{ϵ}_f}}{{ϵ}_f^2}\right)=\int dV{\delta }_{\mathrm{D}}^2{\bar{j}}_{{ϵ}_f},$$

$${\overline{ϵ}}_f={ϵ}_f/{\delta }_{\mathrm{D}}.$$

$${\delta }_{\mathrm{D}}(\theta ,\varphi ,{\Omega }_{b}t,i)=\frac{1}{\Gamma (\theta )\left(1-{\beta }_{\Gamma }(\theta ){\hat{\mathbf{n}}}_v\cdot {\hat{\mathbf{u}}}^{\prime }\right)},$$

$$\mathbf{u}$$

$${\hat{\mathbf{u}}}^{\prime }$$

$$dV\approx \Delta (\theta ,\varphi )dA=a^2\Delta (\theta ,\varphi )R(\theta )\sin \theta \sqrt{R(\theta {)}^2+(dR/d\theta {)}^2}d\theta d\varphi .$$

$$L_{\mathrm{SR}}$$

$$\varsigma (\theta ,\varphi )$$

$${\hat{\mathbf{r}}}^{\prime }={\Lambda }_i{\Lambda }_{{\Omega }_{b}t}\hat{\mathbf{r}},$$

$${\Lambda }_{{\Omega }_{b}t}$$

$$\text{MSP (primary): }{\mathbf{r}}_{\mathrm{NS}}(t)$$

$$\text{Companion (secondary): }{\mathbf{r}}_{\mathrm{c}}(t)$$

$$r_i(t)=\frac{r_i(1-{ϵ}^2)}{1+ϵ\cos \left({\Omega }_{\mathrm{b}}t\right)},r_i=r_{\mathrm{c}},r_{\mathrm{NS}}$$

$${\mathbf{r}}_{\mathrm{NS}}^{\prime }(t)$$

Full text

radiation mechanisms: non-thermal — pulsars: individual (J1023+0038, B1957+20) — binaries: eclipsing — X-rays: binaries

1 INTRODUCTION

This small subset of radio and $\gamma$-ray MSP binaries in tight circular orbits with low-mass companions are useful astrophysical laboratories for the physics of pulsar winds and relativistic shock acceleration. They are relevant to not only striped pulsar winds but also to the physics of Poynting-flux-dominated relativistic outflows in active galactic nuclei and gamma-ray bursts. Observer-dependent high-energy light curves and spectra advance constraints on the underpinning physical phenomena due to observer-dependence of sampling, via orbital modulations, of the emission region in a viewing geometry constrained by radio and optical determinations of the binary mass functions. Prior to the launch of the Fermi Large Area Telescope (LAT; Atwood et al., 2009) $\gamma$-ray observatory, only three such binaries were known in the Galactic Field. The first of these is the original “black widow” B1957+20 (Fruchter et al., 1988), a $1.6$ ms MSP orbited by a stellar companion of mass $M_{\rm c}\gtrsim 0.02M_{\odot}$ with a binary period of $9.17$ hours. It is now well-established that old “recycled” MSPs emit $\gamma$-rays up to several GeV with a spectrum similar to many young pulsars (Abdo et al., 2013) and that prolific e+e- pair cascades (Sturrock, 1971) must occur in the MSP magnetospheres (Venter et al., 2009); some of these pairs are advected into the relativistic pulsar wind that then interacts with the companion star and wind. Follow-up observations, predominantly in the radio band, of Fermi LAT unidentified sources have expanded the binary population to over 30111https://confluence.slac.stanford.edu/display/GLAMCOG/Public+List+of+LAT-Detected+Gamma-Ray+Pulsars in the Field, bringing their total number to over 70 known when including those residing in globular clusters (Manchester et al., 2005) .

Precision radio timing of the binary MSPs accounting for orbital Doppler wobbles yields the pulsar binary mass function andan estimate of the minimum mass of the companion and semi-major axis of the orbit, typically $\sim 10^{11}$ cm. Empirically, the known population of rotation-powered MSP binaries are loosely segregated based on the minimum companion mass $M_{\rm c}$ (Roberts, 2011): black widows (BWs) with minimum companion masses $M_{\rm c}\lesssim 0.05M_{\odot}$ that may be degenerate, and the rarer redbacks (RBs) with non-degenerate companions $M_{\rm c}\gtrsim 0.1M_{\odot}$. These MSP binaries, colloquially termed “spider” binaries, are ancient with characteristic ages $>$Gyr. They “devour” and destroy their low-mass companion by accretion followed by ablation and mass loss exacerbated by the pulsar wind. Recently, some RBs have been observed to transition between a rotation-powered pulsar state and a low-mass X-ray binary accretion state (Archibald et al., 2009; Papitto et al., 2013), confirming the association between these source classes. The behavior of these transitional objects is complex, exhibiting poorly understood transitions in X-ray luminosity and accretion states (Linares, 2014; Bogdanov et al., 2015). The focus of this paper is modeling the conceptually clearer rotation-powered state of BWs and RBs where the total energy budget is constrained by the pulsar rotational energy loss rate $\dot{E}_{\rm SD}$.

In the standard narrative envisioned for BWs and RBs, an intense $\dot{E}_{\rm SD}\sim 10^{34}-10^{35}$ erg s*-1* MSP pair plasma “striped” wind reprocesses in an intrabinary shock that irradiates the tidally-locked companion, preferentially heating the facing side. The companion exhibits orbitally-modulated optical light curves interpreted as a convolution of ellipsoidal variations due to the companion being nearly Roche lobe-filling and anisotropic photospheric emission due to pulsar irradiation. Besides heating by irradiation, particle acceleration beyond TeV energies can take place in the relativistic magnetized intrabinary shock (Harding & Gaisser, 1990; Arons & Tavani, 1993), whose pressure support derives from either the companion wind or a magnetosphere. The companion wind matter also generates radio eclipses of the MSP at orbital phases where the companion is between the pulsar and observer. In this picture, eclipses in the radio pulsations of the MSP in some systems can assist in constraining the shock and wind geometry as we demonstrate in this paper. A photometric light curve model of the orbital optical variations of the companion can be used to constrain the system inclination and irradiation efficiency (e.g., Breton et al., 2013). If radial velocities of the companion can be measured spectroscopically, the companion mass function constrains orbital parameters in the usual way. Combining the radio pulsar mass function, optical companion mass function and modeled system inclination then yields the complete orbital solution of the system, including the neutron star mass. This procedure has been applied to a handful of systems yielding heavy neutron star masses $M_{\rm MSP}$ well-above the canonical $1.4M_{\odot}$ (e.g., van Kerkwijk et al., 2011). Such massive stars constrain theories of the nuclear equation of state, marking these systems as attractive targets for the Neutron star Interior Composition Explorer (NICER, Arzoumanian et al., 2014) due for launch in 2017 with an energy range of $0.2$ – $12$ keV and sensitivity about twice that of XMM-Newton.

The scrutiny of BWs and RBs can help uncover the largely unknown physics of pulsar winds in MSPs, and aid in understanding where the transition occurs from a magnetic flux-dominated to a particle-dominated flow. Unlike the Crab Nebula whose termination shock or inner knot is $\sim 10^{15}-10^{17}$ cm away from the pulsar, the companions in BWs provide a fixed target at a distance only $\sim 10^{11}$ cm from the pulsar with much higher magnetic fields realized than in PWNe shocks. Kinetic-scale magnetic dissipation (i.e., shock-driven reconnection, e.g., Sironi & Spitkovsky, 2011a, b; Lyutikov et al., 2016) in the shocked pulsar wind is a probable acceleration process for leptons if the pulsar wind magnetization $\sigma$, the ratio of magnetic to pair plasma particle kinetic energy density, is larger than unity – however too large a $\sigma$ may preclude the existence of the observed shock. Conversely, if $\sigma$ is small, a more conventional diffusive shock acceleration (DSA) may be the energization mechanism but is likely less efficient due to the oblique shock geometry. Moreover, leptons also may or may not be accelerated in the far upstream pulsar wind, although this scenario is under contention (Lyubarsky & Kirk, 2001; Aharonian et al., 2012; Zrake, 2016). There might be feedback between the intrabinary shock and the upstream wind content, as well (Derishev & Aharonian, 2012).

Unlike massive TeV binaries such as B1259–63, the intrabinary shock in some BWs envelopes the companion rather than the pulsar since the pulsar wind ram pressure dominates. Although many physical processes employed in BW and RB models mirror those invoked in massive TeV binaries (Tavani & Arons, 1997; Dubus, 2013) the shock and orbital geometries are qualitatively different as depicted in Figure 1 . However, as we argue in this paper, many RBs and transitional systems support the interpretation of being “inverted”, where the interaction shock orientation is reversed. It then bows around the pulsar outside the light cylinder of radius $R_{\rm LC}=cP_{\rm MSP}/(2\pi)$, where $P_{\rm MSP}$ is the MSP spin period, rather than around the companion. Such inversions can be envisaged as a state preceding or following accretion in transitional systems, . The shock geometry may significantly impact models of orbitally-modulated high-energy emission as well as the shrouding of the MSP in radio.

In this paper, we construct semi-analytical geometric models for radio eclipses and the Doppler-boosted orbitally-modulated X-ray light curves to constrain the geometry and orientation of the intrabinary shock. Previous analyses of MSP “spider” binaries have largely focused on BW B1957+20, and principally its radio aspects (e.g., Rasio et al., 1989; Tavani & Brookshaw, 1991), leaving the double-peaked (DP) X-ray orbital modulation found in many BWs and RBs (see §2.1 and Table 2.1) unmodeled. We note the parallel and independent work by Romani & Sanchez (2016) for the DP modulation shares some conclusions of this work. We focus on BW B1957+20 and RB J1023+0038 as representative systems, but our framework is generically applicable to other MSP binaries, setting the stage for a future population analysis as well as aiding future models of particle transport in the shock. Such advances will go beyond previous analyses invoking inverse Compton by Bednarek (2014) and Wu et al. (2012), and aid target selection for orbitally-modulated high-energy emission for Fermi LAT and the planned Čerenkov Telescope Array (CTA). In §2 we present interpretations of recent observational developments and construct machinery to constrain the intrabinary shock with a simple geometric model for frequency-dependent radio eclipses. We explore the implications of recent observations and results for shock mixing in §3. In §4 we develop a semi-analytical model for the Doppler-boosted orbitally-modulated X-ray emission, developing synthetic light curves that will be useful for a planned model fitting study.

2 CONSTRAINING THE INTRABINARY SHOCK IN BLACK WIDOW AND REDBACK SYSTEMS

2.1 Overview and Idealizations

As is generically the case for collisionless astrophysical shocks, there are two components separated by the contact discontinuity for the intrabinary shock in MSP binaries like B1957+20 and J1023+0038: a relativistic pair shocked pulsar wind and an ionized shocked companion component (Phinney et al., 1988; Eichler & Usov, 1993). Such a gross structure must exist and is borne out in hydrodynamic (Bucciantini, 2002; van der Swaluw et al., 2003; Bosch-Ramon et al., 2012; Dubus et al., 2015) and relativistic magnetohydrodynamic (RMHD, e.g., Bucciantini et al., 2005) simulations of pulsar wind shocks in various contexts. Relativistic plasma/magnetic turbulence and electron acceleration, likely mediated by magnetic reconnection, DSA, or energization mediated by shear flows (Liang et al., 2013) if mixing between components is low, will occur near the site of the shock contact discontinuity, leading to high-energy emission. For DSA, however, it is widely accepted that oblique relativistic magnetized shocks are less efficient accelerators than parallel shocks and may lead to spatial dependence in the acceleration in the bowed “head” of the intrabinary shock. Due to the disparate length scales of the shock and gyroscale acceleration, such particle acceleration by current kinetic-scale simulations (e.g., particle-in-cell codes) cannot be computed in a self-consistent manner over the large length scales of the shock. Developing an expedient formulation to empirically diagnose the spatial character of such acceleration from high-energy spectroscopically phase-resolved light curves is a planned goal for future work.

To furnish insight to the reader and underpin the framework developed in this paper that will be utilized in the future, we begin by briefly discussing select multiwavelength phenomenology that constrain the shock, physics and geometry in BWs and RBs, with a focus on B1957+20 and J1023+0038. Such an assessment of extant empirical conclusions is critical for the development of a self-consistent, unified model of BWs and RBs as we attempt in this paper.

Phase zero of the orbit in most observational contexts for MSP binaries is defined where the pulsar passes the ascending node, with orbital phases $0.25$ and $0.75$ the superior conjunction (SC) and inferior conjunction (IC) of the MSP, respectively — in this article, phase zero is defined as SC where the MSP is behind the companion for the observer as this is a natural choice for radio eclipses. For some MSP binaries where the pulsar is totally shrouded in the radio and no radio ephemeris is available, the IC phase is associated with the global optical maximum of the irradiated stellar companion such as for J2339.6–0532 (Romani, 2015; Salvetti et al., 2015).

2.1.1 X-ray Phenomenology and Interpretation

Many rotation-powered BWs and RBs show evidence for orbitally-modulated X-ray emission, likely due to synchrotron cooling of relativistic electrons and positrons at an intrabinary shock in a turbulent and relatively high $\lesssim 50$ G magnetic field anticipated just upstream of the shocked pulsar wind (see Eq. [1]). The BW or RB systems, without detectable accretion or disks in the rotation-powered state, have an inherently different origin of X-ray emission than for dipping/eclipsing LMXBs (e.g., Parmar et al., 1986) where accretion power may dominate. The emission typically has a strong non-thermal power-law component, with relatively flat photon indices $\Gamma_{X}\sim 1-1.5$ (e.g., Roberts et al., 2015) that implies relatively hard underlying electron distributions $p\approx 2\Gamma_{X}-1\lesssim 2$ and efficient acceleration.

The significant orbital modulation in spiders is generally not strongly energy-dependent (e.g., Bogdanov et al., 2011, 2014b; Tendulkar et al., 2014) in phase-resolved spectroscopic studies in the soft X-ray band, thus photoelectric absorption is disfavored as the principal modulating mechanism in most sources (J2339.6–0532 may be an exception, cf.  Yatsu et al., 2015). Deep soft X-ray observations of a large subset ($\sim 25\%$) of MSP binaries in the rotation-powered state exhibit strong DP orbitally-modulated X-ray fluxes, most often centered at IC, with BW B1957+20 unique by being centered at SC (cf. Table 2.1). Some other BWs with shallower observations also show hints in photon counts for SC-centered DP emission (e.g. J1810+1744, Gentile et al., 2014), suggesting that the preponderance of IC-centered DP emission may be an observational selection bias for brighter RBs and transitional systems. Observations by NuSTAR (Harrison et al., 2013) of J1023+0038 in its rotation-powered state also reveal the nonthermal orbital modulation extends to at least $79$ keV (Li et al., 2014), with hints for IC-centered DP emission even in the $3$–$79$ keV harder X-ray band (Tendulkar et al., 2014). Such hard X-ray phenomenology is also present in other spiders (M. Roberts, private communication).

Simple occultation by the companion as invoked for J1023+0038 by Bogdanov et al. (2011) cannot naturally explain the DP light curve structure centered around IC as observed by Archibald et al. (2010) and Tendulkar et al. (2014), Moreover, many of these systems have binary inclinations well away from edge-on ($34^{\circ}$ to $53^{\circ}$ for J1023+0038, Archibald et al., 2009), requiring the emission region (and intrabinary shock) relatively close to the companion in the occultation model. Although the X-ray emitting and radio eclipse regions need not be coincident, the large $>50\%$ orbital fraction of radio shrouding of the MSP suggests the plasma is not well-confined near the companion. However it is difficult to envision a plausible and relatively stable hydrodynamic scenario where a shock exists near the companion $L_{1}$ point but other plasma is shrouding the pulsar $\gtrsim 50\%$ of the orbit but generally not at pulsar IC for such low $\lesssim 55^{\circ}$ inclinations. Moreover, for an X-ray emission region close to the companion, occlusion also innately leads to a DP structure that has a peak separation of $\sim 0.5$ that is too wide for any observed BW or RB. Unlike eccentric TeV binaries, the DP light curves in circularized BWs and RBs also cannot be explained by dynamical changes of shock radius and particle cooling between periastron and apastron (e.g., Tavani & Arons, 1997).

We argue in this paper that geometric Doppler boosting of emission along an intrabinary shock, either bowed toward or away from the companion, can naturally explain the DP light curve structure centered at SC or IC, respectively. Then, the phase centering of the DP structure is a key discriminant of the shock orientation and system state. In addition, the light curve structure serves as a probe of shock geometry, particle acceleration, and shock mixing. The bulk Lorentz factor that predicates the Doppler boosting is critically dependent on the level of mixing between the relativistic e+e- wind and the shock-heated ionized companion matter, that is, the baryon loading of the flow. For a striped wind of magnetization $\sigma$ where the shock approximately lies around the line joining the two stars, the striped pulsar wind field orientation relative to the shock normal is critical for particle acceleration (e.g., Sironi & Spitkovsky, 2011b; Summerlin & Baring, 2012). For a striped wind that is envisioned as parallel slabs of alternating field orientation, the shock geometry is quasi-perpendicular at the nose with the highest compression ratio, transforming smoothly to quasi-parallel at the flanges with a lower compression ratio. This spatial dependence of the compression ratio, relativistic shock obliquity, along with higher particle resident time near the stagnation point (the fast cooling locale in Figure 1), should inherently influence the local particle acceleration, cooling and emergent radiation depending on what shock locales the observer line-of-sight samples as a function of orbital phase. However, a detailed exploration in a self-consistent geometry with a transport model for leptons along the shock is deferred to a future paper. For our present study, we focus on the gross DP structure of the light curves in different geometries that can easily be adapted for different sources and energies.

It can be shown that the equatorial upstream wind magnetic field magnitude $B_{\rm w}$, dominated by the toroidal component at large cylindrical radii $r_{\rm s}\gg R_{\rm LC}$ from the pulsar, is

[TABLE]

For a well-defined MHD shock to develop, the magnetization must attain $\sigma\ll 1$ upstream of the shock, either by shock-mediated reconnection (Sironi & Spitkovsky, 2011b) very near the shock precursor, or other kinetic-scale dissipation processes far upstream. Neglecting any baryonic mass loading, the condition $\sigma=B^{2}/(4\pi\,n_{\rm e,MSP}\langle\gamma_{\rm w}\rangle m_{e}c^{2})\lesssim 1$ with Eqs. (1)–(LABEL:neMSP) implies a mean Lorentz factor $\langle\gamma_{\rm w}\rangle$ for an isotropic pair wind,

[TABLE]

where ${\cal M}_{\pm}\gg 1$ is the pulsar pair multiplicity of the Goldreich-Julian rate, cf. Eq. (LABEL:GJchargerate). Following attaining $\sigma\ll 1$, the magnetic field in the shocked pulsar wind field $B_{\rm s}$ then scales as $B_{\rm s}\sim 3\sqrt{\sigma}B_{\rm w}$ in the ultrarelativistic perpendicular shock limit (Kennel & Coroniti, 1984). However, the magnetic dissipation processes upstream may convert or destroy the striped wind morphology such that the shock may be quasi-parallel in the proper frame. A containment argument, based on the observed X-ray power law provides a rudimentary lower bound on $B_{\rm s}$ – the Larmor radius $r_{\rm L}$ of electrons in the shock must be smaller than about $1\%$ of the orbital length scale $r_{\rm L}\lesssim 0.01a\sim 10^{9}$ cm. Then, assuming emission at the critical synchrotron dimensionless energy $\epsilon_{c}=3B_{s}/(2B_{\rm cr})\gamma_{\rm e}^{2}$ with $B_{\rm cr}\approx 4.414\times 10^{13}$ G and electron Lorentz factor $\gamma_{\rm e}$, for an observed power law extending to energy $\epsilon_{X,\rm max}$ in units of $m_{e}c^{2}$,

[TABLE]

where we have neglected factors of roughly unity associated with Doppler shift of energies corresponding to mildly relativistic bulk speeds along the shock. Therefore power laws extending up to $\epsilon_{X,\rm max}\approx 0.15\approx 80$ keV$/(m_{e}c^{2})$ observed by NuSTAR for J1023+0038 advance $2\lesssim B_{\rm s}\lesssim B_{\rm w}\sim 200$ G in the relativistic magnetized shock if $r_{s}\sim 10^{10}$ cm, which implies radiating electron Lorentz factors of order $10^{5}$–$10^{6}$ For other spiders where observations at energies above the classical soft X-ray band are not available, the field magnitude is still greater than about one Gauss, orders of magnitude larger than those in plerions. We consider implications of these bounds on the shock in §3.

2.1.2 Radio Phenomenology

Orbital eclipses of the MSP’s radio pulsations are a common feature in many BWs and RBs in the rotation-powered state. Observed orbital eclipse fractions $f_{E}$ are ordinarily $f_{E}\sim 5-15\%$ for BWs, and typically much larger for RBs, increasing in low radio frequency bands. For example, PSR J1023+0038 eclipses for less than $5\%$ at $3$ GHz to over $\sim 60\%$ of an orbit at 150 MHz (Archibald et al., 2009, 2013). Some BWs also have extensive eclipses. There appears to be a dichotomy in the relative stability of eclipses – for some BWs like B1957+20 eclipses near SC are generally stable orbit-to-orbit, while sporadic mini-eclipses are seen in some other systems particularly those systems with larger eclipse fractions (e.g., Archibald et al., 2009; Deneva et al., 2016). However even in these erratic systems with mini-eclipses, the pulsar is generally unshrouded at IC in relevant bands. A standard decomposition of $f_{E}$ into symmetric and antisymmetric parts about SC is attainable as a function of observer frequency $\nu$. Frequency dependence of the eclipse fraction asymmetry is standard, with larger asymmetry in ingress-egress delays at lower observing frequencies, e.g. PSR B1957+20 (Ryba & Taylor, 1991; Stappers et al., 2001) and J1023+0038 (Archibald et al., 2009, 2013). At the highest radio frequencies $\nu$, the antisymmetric part of $f_{E}$ is typically small compared to the symmetric part.

For B1957+20 and other systems, the symmetric part of these eclipses encompass inferred length scales that are significantly larger than $R_{*}$ for a fully Roche lobe-filled companion, even for $\sin i\approx 1$. No eclipses by the companion are expected if $i<90^{\circ}-\arcsin(R_{*}/a)$, but many systems with eclipses have well-constrained inclinations and companion sizes which violate this inequality. Therefore, eclipses must be predicated on plasma within the system and/or a secondary magnetosphere. Eclipses typically exhibit large plasma dispersion measures before the coherence in the timing solution of pulsations is lost, likely due to absorption rather than scattering (Roy et al., 2015); continuum eclipses of the pulsar are also seen in some systems at low frequencies e.g., for BW B1957+20 (Fruchter & Goss, 1992) and RB J2215+5135 (Broderick et al., 2016) with a scaling $f_{E}\propto\nu^{-0.4}$.

There are a panoply of potential eclipse mechanisms (cf. Michel, 1989; Eichler, 1991; Gedalin & Eichler, 1993; Thompson et al., 1994) depending on physical parameters realized in the intervening plasma. Cyclotron absorption has been posited in B1957+20 (Khechinashvili et al., 2000) but relatively little Faraday rotation is seen, consistent with a $1-10$ G mean magnetic field magnitude in the eclipsing medium (Fruchter et al., 1990), not inconsistent with Eq. (3) since the eclipsing medium consists of the ionized companion wind as well. Moreover, it is now known that the companion in B1957+20 is likely non-degenerate (Reynolds et al., 2007). Excess delays, consistent with plasma dispersion, generally show that the average free electron column density rises sharply from $\langle n_{e}\rangle d\sim 10^{15}$ cm*-2* to $10^{18}$ cm*-2* at phases deep into the eclipse (Ryba & Taylor, 1991; Stappers et al., 2001) for BWs, for $d\sim a$ the line-of-sight column depth, but it is anticipated that there is also clumping near the shock contact discontinuity. This $\langle n_{e}\rangle$ is much higher than implied by Eq. (LABEL:neMSP), therefore the companion wind must have some influence. Whatever the mechanisms for eclipses, the momentum flux balance between the pulsar wind and a companion wind or magnetosphere defines a geometric volume of plasma through which the MSP is eclipsed, bounded by the shock surface (gray curves depicted in Figure 1).

Consequently, we advance that the dichotomy of eclipse phenomenology is the orientation of the shock surface germane to the X-ray light-curve phasing in Table 2.1. For the SC-centered DP phase centering where the shock is bowed around the companion, as for BW 1957+20, the relative stability and small $f_{E}$ are consistent with this picture. Contrastingly, for IC-centered X-ray phasing where the shock is orientated around the pulsar, larger and more erratic eclipses are expected where the companion wind can enshroud the pulsar, and is necessarily turbulent for the obligatory angular momentum loss. The radio optical depth, as well as the shock orientation depend on the companion wind mass loss rate. This can be very low or substantial through evaporation or quasi-Roche lobe overflow (e.g., Bellm et al., 2016), respectively, but is poorly understood.

For the IC-centered scenario, canonical Roche lobe overflow at the characteristic ion sound speed cannot be a wind source since the circularization radius $R_{\rm circ}$ must be larger than the shock radius $R_{0}$ (measured from the MSP), or the system will be predisposed to a disk-state (Frank et al., 2002). Moreover, for the radio pulsar state, $R_{0}$ must exceed the light cylinder scale, that is, $R_{0}>\mbox{Max}(R_{\rm circ},R_{\rm LC})$. This then favors an evaporatively-driven quasi-Roche lobe overflow supersonic wind model for rotation-powered states. The mass loss must be low enough to escape IR/optical detection. The scenario is somewhat fine-tuned such that the companion wind is fast enough to inhibit a disk, while dense enough such that angular momentum losses owing to turbulence are sufficient for gravitational influences to overpower the pulsar wind. Such turbulence may also be driven by the radio absorption that predicates the eclipsing mechanism. Accordingly, $R_{0}/{\cal{R}}\sim(\dot{m}_{g}c^{2}/\dot{E}_{\rm SD})^{2}\gg 1$ where ${\cal{R}}=2GM_{\rm MSP}/c^{2}$ is the Schwarzschild radius of the MSP and $\dot{m}_{g}$ is the gravitationally-captured wind’s mass rate. However there are stability concerns – the pulsar termination shock that arrests accretion flow and shrouds the pulsar and delineates the eclipsing medium may only be pushed out to a modest $2-100$ multiples of pulsar light cylinder radius $R_{\rm LC}\ll a$ (Ekşİ & Alpar, 2005; Linares, 2014) unless a feedback mechanism is operating. We show in an upcoming paper that if there is a feedback mechanism operating in RBs for the mass-loss rate from the companion, then such autoregulation may permit the shock to be stable much farther from the light cylinder out to orbital length scales for rotation-powered disk-free states (Wadiasingh 2017, in prep). A detailed discussion of the poorly-understood nuances of the irradiated companion mass loss, stability, and eclipse mechanisms is beyond the scope of this paper.

2.2 Geometric Constraints by Radio Eclipses

Here we explore what constraints on the shock parameters can be gleaned from just the geometry of eclipses as a function of inclination. This requires an a priori model for the shock geometry. For cases where the shock is orientated around the companion, two formalisms have been invoked for the free electron density underpinning the eclipses: an optically thin low-density plasma tail from the companion that spans several orbital semi-major axis length scales $a$, and an optically-thick model of much higher local free electron density in the intrabinary shock and shocked companion wind (Rasio et al., 1989, 1991; Tavani & Brookshaw, 1991). We adopt the optically-thick formalism, as this leads to constraints that are upper limits on $R_{0}$ for a given $i$ which we apply to B1957+20 as a test case in § 2.2.1. That is, for a model shock surface geometry, the transition from optically thick to thin may be parameterized in terms $\theta_{\rm max,R}$ for a given radio band. This is conceptually similar to radius-to-frequency mapping used in the study of pulsed emission in radio pulsars (Komesaroff, 1970; Cordes, 1978). Small values of $\theta_{\rm max,R}$, corresponding to high radio frequencies, sample regions of the shock closer to the shock head. Asymmetry of eclipses is interpreted as Coriolis influences on the shock, skewing it by an angle or sweeping-back a cometary tail; the latter is explored in the Appendices.

At this stage, we do not attempt to self-consistently model the shock geometry and parameters from, for instance, generic covariant MHD jump conditions (Double et al., 2004).

The general problem of an arbitrarily-shaped region occulting a source involves computational geometry techniques that may require inefficient ray casting or tessellating grids. To make the problem more analytically expedient, we assume an azimuthally symmetric form for the intrabinary shock with radial function $\{R(\theta)\,|\,\theta\in(0,\theta_{\infty})\}$ along the axis of symmetry of the bow-shaped shock, with $R(\theta=0)=R_{0}$, and generalize this approach to approximate the swept-back tails due to orbital motion. Azimuthal symmetry of the shock is expected to be an acceptable approximation for the intrabinary shock locale in the vicinity of the shock stagnation point, although overall the shock angle may be skewed relative to the line joining the two stars yielding the extended-egress delayed-ingress radio eclipse phenomenology. This approximation to the shock structure should especially be good when restricted to the ingress portion of eclipses about SC that are sharp and well-defined, unlike those that are more diffuse at egress and that are likely contaminated by plasma from the cometary-tail. In Appendix A, we develop an analytical formalism for eclipses of a point source by an arbitrary azimuthally-symmetric surface orbiting the source around a common barycenter, as well as parameterize the asymmetry of eclipses due to the swept-back tail for finite flow velocities and wind accelerations.

In the highly-radiative supersonic limit, azimuthally-symmetric analytic purely-hydrodynamic bow shock forms that assume nonrelativistic, momentum-dominated winds neglecting gravity, through the balance of ram pressures, can be found in Wilkin (1996) and Canto et al. (1996). Here the companion shock, contact discontinuity, and shocked/deflected pulsar wind are roughly spatially coincident compared with the length scale $a$, although not necessarily well-mixed, and internal pressure contributions to the momentum-flux tensor are neglected. If the flow is not supersonic, the analytic forms are significantly narrower than simulations of hydrodynamic shocks with finite Mach number (e.g., van der Swaluw et al., 2003). Relaxing the highly-radiative limit or introducing mass loading, parameterizations for azimuthally-symmetric shocks can also be found in Gayley (2009) or Morlino et al. (2015), and generically result in increasing the shock opening angle.

Since the physics and geometry of the pulsar wind and induced companion wind or magnetosphere are largely unknown, we utilize a generic two-isotropic-colliding-winds analytic solution Canto et al. (1996) for the intrabinary shock geometry that surveys varied shock asymptotic angles through a simple parameterization. This is readily apparent from the gray contours in Figure 1 and qualitatively resembles hydrodynamic simulations of Tavani & Brookshaw (1991) for B1957+20. As explored in §2.2.2, the geometric form is immaterial for radio eclipses for the case where the shock is orientated around the pulsar. However, for eclipses where the shock is around the companion, the prescribed geometry is consequential for the parameter constraints. For this reason, we consider an alternative parallel-wind geometry of Wilkin (1996), a standard bow shock, in Appendix LABEL:parawindappendix which has significantly narrower asymptotic shock angles and is the $R_{0}\ll 1$ limit of the two-wind solution near the shock head. We limit ourselves to axisymmetric forms for simplicity – prescriptions for the shock geometry that include nonaxisymmetric distortions by Coriolis effects are found in Parkin & Pittard (2008).

The radial function for the colliding isotropic winds parametrizing the shock geometry is

[TABLE]

and $\theta\in(0,\theta_{\infty})$ the polar angle defining the shock, with zero taken as the line separating the two stars and $\theta_{\infty}$ the asymptotic shock angle,

[TABLE]

where $\eta_{\rm w}$ is the ratio of the two wind ram pressures and $\theta_{1}$ is an implicitly defined function of $\theta$ and $\eta_{\rm w}$. Explicitly, $R_{0}$ is related to $\eta_{\rm w}$ by

[TABLE]

We caution that the physical interpretation of $\eta_{\rm w}$ may be misleading, especially for scenarios where the shock wraps around the pulsar where gravitational influences and unknown wind anisotropies are salient. Otherwise, when the shock is orientated around the companion, as for B1957+20 in §2.2.1, $\eta_{\rm w}$ may be connected to the pulsar’s parameters if the wind of the companion is induced, e.g., Eqs. (10)–(11) in Harding & Gaisser (1990) and Eq. (7) in the next Section.

In §4, we only employ the geometry of these analytical shocks rather than their physical velocity and density profiles, since the companion and pulsar shock components may not be well-mixed. Generically for such a pressure-confined flow, surface mass density for the bow shocks is highest near the stagnation point and slowly decreases farther down, while tangential velocity increases approximately linearly with the shock polar angle parameter $\theta$. These are readily apparent from a Taylor series expansion of analytical expressions in the thin-shell limit. At the stagnation point, the tangential velocity approaches a small value, if internal pressure contributions are small – this is indeed the behavior observed in nonrelativistic hydrodynamic simulations (e.g., Figure 3 of Reitberger et al., 2014). The tangential velocity $v_{\rm s}$ and surface density $\Sigma_{\rm e}$ variation with $\theta$ for hydrodynamic bow shocks can be shown to follow $v_{\rm s}\propto\theta$ and $\Sigma_{\rm e}=\Sigma_{0}(1+w\theta^{2})$ for $\theta\lesssim 1$ with $w<0$, where $|w|\ll 1$ is a constant, but the exact normalization of these generic forms depends critically on the wind pressures and velocities that are unknown. This generic form, generalized to relativistic velocities, is utilized in §4 for the semi-analytic DP light curve synthesis. However, the hydrodynamic density profile is probably an inaccurate proxy for the spatial distribution of particle acceleration and cooling; thus a more general procedure is also described in §LABEL:partaccdiag.

2.2.1 Application to PSR B1957+20, an SC-centered Spider

Using the method developed in Appendix A, we compute the axisymmetric eclipse fraction $f_{E}$ for the shock geometry in Eq. (4) using Eq. (A21). This is conditional on the crucial parameter $\theta_{\rm max,R}$, where the medium transitions from optically thick-to-thin for a given observing frequency. Note that $R_{0}$ and $i$ are the geometric parameters that are independent of observer frequency, therefore $\theta_{\rm max,R}$ is the only parameter in the model that connects to the frequency dependence of symmetric eclipses.

In Figure 2, we display the eclipse fraction $f_{E}$ dependence on $R_{0}$ of PSR B1957+20 with a fixed mass ratio $q=69.2$ for various inclination angles consistent with companion light curve models (Reynolds et al., 2007; van Kerkwijk et al., 2011) as well as $\zeta\approx i\approx 85^{\circ}$ found from Johnson et al. (2014). The axisymmetric computations for constraining $R_{0}$ should be more accurate for the eclipses at $\nu\gtrsim 600$ MHz that are largely symmetric about SC; this is indicated by the red region in the panels. We neglect the minor degenerate observational coupling between mass ratio and inclination angle due to uncertainties in the optical mass function. The four panels depict successively larger values of the parameter $\theta_{\rm max,R}$ left to right. In the leftmost two panels the value of $\theta_{\rm max,R}=\{\pi/3,\pi/2\}$ is independent of $R_{0}$, while the rightmost two panels impose values that depend on $R_{0}$, e.g. through Eq. (5). In particular, the prescription $R(\theta)\sin\theta=4R_{0}$ selects the value of $\theta_{\rm max,R}$ for a given $R_{0}$ such that the transverse shock length scale is $4R_{0}$. On such a scale far from the shock head, hydrodynamic instabilities will likely develop that may break the assumptions in our rudimentary model. The rightmost panel chooses an exceptionally large value of $\theta_{\rm max,R}=0.9\theta_{\infty}$ that serves as an extreme limit for what $\theta_{\rm max,R}$ may be and represents a very substantial occluding volume. Such large values of $\theta_{\rm max,R}>\pi/2$ include regions well beyond the head of the shock, and should not be axisymmetric due to Coriolis effects and instabilities. Despite that well-defined and sharp symmetric eclipses are not expected from such large values of $\theta_{\rm max,R}$, these values are included for completeness. Similarly, values much smaller than $\theta_{\rm max,R}<\pi/3$ require rather flat shocks for a fixed $f_{E}$ tending to $R_{0}\rightarrow 0.5$ well past the $L_{1}$ point, also a rather unreasonable scenario that requires $\eta_{\rm w}\sim 1$ and is in tension with the X-ray DP light curve peak separation in §4.

For fixed inclination and $f_{E}$, it is clear from Figure 2 that larger values of $\theta_{\rm max,R}$ are compatible with smaller values of $R_{0}$. For large inclinations near edge-on, there are clearly constraints on $\theta_{\rm max,R}$; in particular for $i=85^{\circ}$, $\theta_{\rm max,R}\lesssim\pi/2$ for $R_{0}>R_{*}$. On the other hand, for $i=65^{\circ}$, the constraint is looser with $\theta_{\rm max,R}>\{\pi/3,\pi/2,R(\theta)\sin\theta=4R_{0},0.9\theta_{\infty}\}$ corresponding to limits $R_{0}\lesssim\{0.31,0.235,0.17,0.14\}$, respectively, for $f_{E}\lesssim 7\%$. The upper limits on $R_{0}$ are modestly more stringent for $\nu>600$ MHz for the same range of $\theta_{\rm max,R}$. Therefore we conclude that for $i=65^{\circ}$, $R_{0}\approx 0.17-0.3$ with a canonical value of $R_{0}\approx 0.235$ for $\theta_{\rm max,R}=\pi/2$ which defines the shock head. This latter value of $R_{0}\approx 0.2$ is employed in §4.2 and seems plausibly compatible with the observed X-ray DP light curve; too large an $R_{0}$ leads to DP light curves that have too-wide of a peak separation as will become apparent in due course. Some geometric realizations for $f_{E}\approx 7\%$, i.e. eclipses at $\nu\approx 600$ MHz, are illustrated in Figure 3 with corresponding numerical values of $\theta_{\rm max,R}=\{\pi/3,\pi/2,R(\theta)\sin\theta=4R_{0},0.9\theta_{\infty}\}$ (columns) and $R_{0}$. Some geometric trends are clearly evident in Figure 3, for instance larger $i$ requiring smaller $R_{0}$ for the same $f_{E}$. Similarly, larger $\theta_{\rm R,max}$ for fixed $i$ allows for lower $R_{0}$.

The upper limits found for $R_{0}$ allow for an estimate of the companion wind pressure if $\dot{E}_{\rm SD}$ is known. We may express $\eta_{\rm w}$, with ${\cal{P}}$ the wind pressure due the companion, as

[TABLE]

where $T_{\rm cold}$ is the isotropic unheated temperature of the stellar companion (i.e., the intrinsic unirradiated radiation pressure from the secondary), $\delta\Omega/(4\pi)$ is the fractional isotropic pulsar wind solid angle subtended by the shock. If a strong magnetosphere is not the principal source of ram pressure against the pulsar wind, then the parameter $\xi$ embodies the fractional energetic efficiency, in units of $\dot{E}_{\rm SD}$, of the induced companion wind generated by unspecified processes. The $R_{0}$ found above are compatible with the thermally-driven wind or companion magnetosphere scenarios of Harding & Gaisser (1990). Numerically, for typical values $R_{0}=0.2-0.3$ consistent with radio eclipses, the ratio of ram pressures $\eta_{\rm w}$ is between $6$ to $18\%$ by inverting Eq. (6). From Eq. (7), we can form an estimate of the energetic efficiency $\xi$ of the induced wind, if the intrinsic or induced magnetic field of the companion is small,

[TABLE]

where the ratio of cold intrinsic stellar to pulsar power can be neglected in BWs and RBs, since it of the order $\sim 10^{-4}-10^{-6}$, and thus $\xi$ is a simple function of stagnation point $R_{0}$. This equation is unphysical for $R_{0}\rightarrow 0.5$ and should not be used in this limit. The solid angle fraction for the canonical shock head $\theta\in(0,\pi/2)$ which may participate in the heating of the companion can be routinely found $\delta\Omega/(4\pi)\approx 3R_{0}^{2}(1+2R_{0})/4$ for $R_{0}\ll 1$ from Eq. (4). Whence, the efficiency of the induced wind $\xi$ is of the order $0.1-5\%$, depending upon the inferred shock standoff $R_{0}$. This efficiency is similar in order-of-magnitude to the hemispherical quiescent induced photospheric heating fraction, i.e. $2\pi R_{*}^{2}\sigma_{B}T_{\rm hot}^{4}/\dot{E}_{\rm SD}\sim 0.4\%$ for B1957+20 where $T_{\rm hot}\approx 8000$ K, as one would expect in a model where the companion’s gas pressure, from a thermally-driven evaporative wind, balances the striped cold pulsar wind.

An order-of-magnitude upper limit to the companion surface magnetic field $B_{*}$ may be found by assuming companion pressure is entirely due to a magnetosphere at the stagnation point (e.g., Eq. (23) of Harding & Gaisser, 1990),

[TABLE]

which is relatively small compared to typical surface fields found in degenerate cores, and comparable to kilogauss fields found in T Tauri stars (Johns-Krull, 2007). In principle, such fields may be detectable in future, but currently requires relatively bright (magnitude $\lesssim 12-14$, C. Johns-Krull, private communication) companions for a sufficient signal-to-noise with high-precision IR/optical spectroscopy of Zeeman broadening on a model atmosphere. Any optical field constraint would be useful in constraining the physics of the pulsar and induced companion winds, as well as to appraise the Applegate (1992) model for a tidally-driven convective dynamo. Indeed, a strong companion magnetosphere would also dramatically alter shock MHD jump conditions, impact the relevant acceleration mechanisms in the shock, and influence particle heating of the companion.

Observe in Figure 3 that the spatial region of the shock the MSP is eclipsed by is markedly different for $55^{\circ}$ and $85^{\circ}$, with the former only sampling regions of the confined flow peripheral to the shock head. This segues to the inclination-dependent numerical computation, using Eq (A21) of the growth rate of the eclipse fraction $f_{E}$ as function of $\theta_{\rm max,R}$, depicted in Figure 4. The appropriate interpretation of Figure 4 should be restricted to the symmetric part of eclipses, i.e. below about $f_{E}\lesssim 0.1$. The curves terminate at $\theta_{\infty}$, and we note the computed $f_{E}<0.5$ is finite and bounded even for the unbound axisymmetric shock geometry, i.e. $\theta_{\rm max,R}\rightarrow\theta_{\infty}$. For high inclinations, where the eclipses largely sample the shock head, the growth rate is approximately linear with $\theta_{\rm max,R}$ contrasting the nonlinear growth rates for lower inclinations which sample the periphery and tail of the shock. If the optical depth $\tau$ due to scattering or absorption by a cross section $\sigma_{\nu}\propto\nu^{-m}$ is given by $1\ll\tau\propto\sigma_{\nu}\Sigma_{\rm e,R}$ for column density $\Sigma_{\rm e,R}$, then the spatial variation of the column density can be probed. Given growth curves $f_{E}\equiv g(\theta)$, for an observed frequency dependence in radio eclipses $f_{E}\sim\nu(\theta)^{-n}$, as observed for B1957+20 and RB J2215+5135 with $f_{E}\propto\nu^{-0.4}$, then $\nu(\theta)\sim g^{-1/n}$ and the spatial distribution of the integrated column for a given optical depth is proportional to $\Sigma_{\rm e,R}(\theta)\sim g(\theta)^{-m/n}$. For instance, if $g(\theta)\sim\theta^{l}$, $n=0.4$ and $m=2$ for free-free absorption in the Rayleigh-Jeans limit, then $d\log\Sigma_{\rm e,R}(\theta)/d\log\theta\sim-5l$, an inclination-dependent line-of-sight column density that can attain a rather steep and sharp profile. This motivates frequency-dependent radio population studies of similarly eclipsing BWs and RBs for inclination-dependent trends of $f_{E}(\nu)$. This is of generally low utility to the X-ray Doppler-boosted emission in §4 since the two populations of electrons are not necessarily concordant, but important for constraining the nature and content of the companion wind by exploring the parameter space of Figure 4.

In Appendix LABEL:parawindappendix, we consider the simpler alternative parallel-wind bow shock geometry of Wilkin (1996) and compute correspondent $f_{E}$ growth curves for B1957+20 and other SC-centered spiders. Such a parallel-wind geometry is self-similar in the $R_{0}\ll 1$ limit of Eq (4) with much smaller shock opening angles. In this geometry $R_{0}$ has no impact on the shock opening angle but only serves as a scaling parameter. Since the companion and pulsar winds are not anticipated to be isotropic but somewhat stronger about the line joining the two stars in the orbital plane, this geometry is a realization where wind anisotropies obligate a much smaller shock opening angle. Such wind anisotropy is expected for an anisotropically-irradiated companion whose evaporatively-driven wind is only influential on the day-side hemisphere of the companion. For this auxiliary geometry, it is found that the constraints on $R_{0}$ are systematically larger than the isotropic-winds geometry, and less sensitive to the value of $\theta_{\rm max,R}$ since $f_{E}$ plateaus in the far-downstream tail of the shock. Since the shock angle is small, an additional upgrade to ballistic tail sweepback and eclipse asymmetry is also explored. This parameterization of eclipse asymmetry in terms of outward wind speed serves as an additional motivation for future radio characterization.

2.2.2 Application to PSR J1023+0038 and Other IC-centered Spiders

For scenarios where the shock surrounds the pulsar and the perpendicular component of the shock $R(\theta)\sin\theta$ is a monotonically rising function of $\theta$, as is the case for all shocks in the limit of small sweepback, the eclipse fraction is independent of shock geometry for any given binary inclination and only depends on the shock’s largest attained polar angle $\theta_{\rm max,R}$. For cases where the shock is swept back due to Coriolis effects which may introduce $\phi$ dependences and break the axisymmetry of the shock, the axisymmetric estimate below is a lower limit on the shrouding fraction (see Appendix A).

If the shock is at an angle $\vartheta_{\rm sb}$ with respect to the orbital angular momentum vector, the eclipse fraction $f_{E}$ is related to $\theta_{\rm max}$ by routine spherical trigonometry,

[TABLE]

We adopt the case $\vartheta_{\rm sb}=\pi/2$ corresponding to the shock axis of symmetry lying in the orbital plane. This yields

[TABLE]

for $\pi/2-i\leq\theta_{\rm max,R}\leq\pi/2+i$ with the lower and upper limits corresponding to $f_{E}=0$ and $1$ respectively. Thus radio eclipses in MSP binaries where X-ray emission is IC-centered constrain the shrouding by the maximum shock polar angle $\theta_{\rm max,R}$ in a model-independent manner, and values of the polar angle are associated with a physical integrated line-of-sight density profile. Note that this formalism allows for axisymmetric shocks that are skewed at an angle $\Delta\phi_{s}$ relative to the line joining the two stars. If $\phi_{\rm in}$ and $\phi_{\rm eg}$ are the orbital phases of ingress and egress eclipses relative to SC, respectively, then the shock asymmetry is directly measurable at a given frequency (corresponding to a $\theta_{\rm max,R}$ in this model) as the antisymmetric part $(1/2)(\phi_{\rm eg}-\phi_{\rm in})\approx\Delta\phi_{s}/(2\pi)$.

The contours corresponding to Eq. (11) are illustrated in Figure 5 with approximate J1023+0038 observational radio eclipse fractions from Archibald et al. (2009, 2013) for a range of inclinations. For parallel-wind shocks, large values of $\theta_{\rm max,R}>\pi/2$ correspond to unphysically long tails $L_{\rm z}\gg 4R_{0}$ where hydrodynamic instabilities are expected to be influential. Moreover, as we show in §4.3, such a geometry is disfavored over an isotropic colliding-winds scenario which generally have larger shock opening angles limited by the asymptotic bow shock opening angle defined by Eq. (5). The maximum values of $R_{0}$ corresponding to $\theta_{\infty}$ are shown in Figure 5, thus, one may constrain the upper limit of the shock opening angle based on the largest stable eclipse fractions observed. For J1023+0038, from eclipses at $150$ MHz of $f_{E,150\rm MHz}\gtrsim 0.6$ this estimate yields, $R_{0}\lesssim 0.4$. We caution that in such IC-centered spiders, the naive physical interpretation of $R_{0}$ constraints in terms of wind ram pressures given by Eqs. (4)–(6) is erroneous due the commanding gravitational influence of the MSP past the $L_{1}$ point, and depends on the specifics of angular momentum loss of the companion baryonic wind.

3 The Downstream Bulk Lorentz Factor $\Gamma$, Shock Mixing, and Baryon Loading

For non-relativistic shocks in Eq. (4), the tangential velocity flow increases approximately linearly with the shock polar angle $\theta$ or symmetry axis $z$ for points close to the stagnation point where the shock is very nearly a spherical cap. Lacking a self-consistent relativistic MHD shock geometry, for simplicity we extend this flow scaling relation to the mildly relativistic regime assuming a scaling for the bulk specific relativistic momentum,

[TABLE]

with $\theta_{\rm max,X}$ corresponding to a characteristic angular scale where the shocked pulsar wind is pressure-confined in this manner

[TABLE]

We adapt this ad hoc prescription Eq. (12)–(13) for the calculation of Doppler factors to synthesize light curves in §4 . This linear spatial dependence of specific momentum is a general characteristic of pressure-confined relativistic flows (e.g., Beskin & Nokhrina, 2006; Komissarov et al., 2007, 2009) . The physical interpretation of the bulk Lorentz factors $\Gamma$ depends on the baryon loading of the shocked pulsar wind by the shocked companion component across the contact discontinuity, and couples to the companion mass loss rate in a nontrivial manner; The relative amount of mixing is also important for the normalization of spectra to the energy budget $\dot{E}_{\rm SD}$ since it augments the leptonic population in the shocked pulsar wind, and is consequential for spectroscopically resolved light curves in §LABEL:partaccdiag but largely unimportant for the geometric considerations in §4 that prescribe values of $\left(\Gamma\beta\right)_{\rm max}$ and other parameters in a fixed energy bin.

We anticipate a partitioned shock morphology in spiders, with a shocked companion wind and pulsar wind, separated by a contact discontinuity. From flux-freezing arguments in an idealized steady-state laminar MHD framework, there is no charged particle transport across the contact discontinuity except perhaps at the stagnation point. This lack of transport seems to be in tension with the assumptions of Romani & Sanchez (2016). RMHD instabilities and kinetic-scale processes alter this overly-simplified picture for mixing. Following Morlino et al. (2015), an estimate of mixing due to kinetic-scale effects may be discerned by consideration of cross-field diffusion of charges. In the Bohm limit of strong turbulence and wave-particle interactions, the perpendicular diffusion coefficient $\kappa_{\perp}$ proceeds at the gyroscale attaining $\kappa_{\perp}\lesssim\kappa_{\parallel}=r_{\rm L}v/3$. Although such strong turbulence is not anticipated in the shocked companion wind, this estimate constitutes a lower limit for the cross-field diffusion timescale $\tau_{{\cal D},\perp}$ of protons,

[TABLE]

where we have taken $v_{\rm p}=10^{9}$ cm s${}^{-1}\gg a/P_{\rm b}$ for the proton speed and assumed $\ell=0.3R_{0}\sim 10^{9}$ cm, the length scale corresponding to the separation between the shock components, roughly that recovered from RMHD simulations (Bucciantini et al., 2005). The value of $\tau_{{\cal D},\perp}$ attained in Eq. (14) is much greater that the typical advection timescale of ions of $10-100$ s in the shocked companion wind, thus mass loading due to diffusion is unlikely to be a significant influence unless $\sigma$ is significantly smaller and particle dominated. Eq. (14) adapted to pairs in the shocked pulsar wind is a factor of $1836$ larger, thus less constraining. The onset of significant baryon loading due to kinetic effects then corresponds to $\langle B_{\rm s}\rangle\lesssim 10^{-3}$ and $\sigma\lesssim 10^{-9}$ in the Kennel & Coroniti (1984) prescription for oblique shocks, values that are largely excluded since they obligate a Lorentz factor of leptons in the shocked pulsar wind $\gamma_{e}\gtrsim 10^{8}$ (which violate the Hillas bound for containment in the shock) to emit at the critical synchrotron frequency corresponding to the $>10$ keV NuSTAR band. The contribution of mixing due to instabilities is difficult to quantify and depend on $\sigma$ in a nontrivial manner, but less germane for $\sigma\ll 1$ as they are likely to be hydrodynamic in nature and occur far from the head of the shock in the tail, at distances scaling with the natural driving scale $R_{0}$.

Additionally, using the magnetic field constraint Eq. (3) and the pair wind density Eq. (LABEL:neMSP), it is simple to form an upper limit for the shocked magnetization parameter $\sigma_{\rm shocked}$ if there is no baryon loading,

[TABLE]

Therefore, modest values of ${\cal M}_{\pm}$ and $\langle\gamma_{e}\rangle$ are sufficient to yield $\sigma_{\rm shocked}\lesssim 1$ which implies that mixing is not compulsory for efficient leptonic acceleration, but probably occurs at some low level. Neutral hydrogen may still cross fields and mix; however, the prolific soft X-ray emission from the shock should ionize any such neutrals close to the shock so that the mass loading length scale is much larger than $R_{0}$. Therefore, the mixing of shock components is insignificant unless small-scale instabilities are prolific near the head of the shock or there are unaccounted for magnetic fields quasi-parallel to the shock normal. Mixing with the companion wind electron/ion plasma also dilutes the energy budget for the shock-accelerated energetic positrons that escape the system. Hence, models that posit MSP binary shocks as significant astrophysical sources of energetic cosmic-ray positrons at Earth assume low mixing (e.g., Venter et al., 2015). Unless the unknown mass loss rate of the companion is small, low levels of mixing are also necessary for the mildly-relativistic bulk Lorentz factors ascribed by Eq. (12)–(13) and in §4 . An observational constraint of mixing on the shocked pulsar wind and underlying leptonic particle distribution also follows from the method outlined in §LABEL:partaccdiag.

4 Geometric Doppler-boosted Orbitally Modulated Synchrotron Emission

4.1 Formalism

We assume local quasi-isotropy in pitch angle of the differential electron number density distribution along the shock in the comoving frame of the bulk flow in a thin emitting region such that $\bar{n}_{e}(\boldsymbol{\bar{r}},\gamma_{\rm e})\rightarrow\Sigma_{\rm e}(\theta,\phi,\gamma_{\rm e})/\Delta(\theta,\phi)$ at each lab-frame coordinate pair $\{\theta,\phi\}$ on the shock, where $\Delta(\theta,\phi)$ is the shock thickness and $\Sigma_{\rm e}$ the differential electron surface density distribution. In reality, the local electron distribution accelerated in relativistic shocks by DSA or prolific magnetic reconnection is expected to be highly anisotropic relative to the mean field direction in the comoving frame. However, for locales near the stagnation point where the residence timescale is large compared to the adiabatic and convection timescales, the flow is expected to be well-isotropized with the field randomly oriented, especially if magnetic dissipation is the principal particle acceleration mechanism. Moreover, for locales where the convection speed is large, Doppler boosting transforms any emission small-scale anisotropy into a narrow cone along the shock, the geometry of which dominates in the observer frame over any smaller-scale anisotropy within the cone due to the underlying electron pitch angle distribution.

The differential volumetric emissivity, in the comoving frame (denoted by bars over variables), is defined by

[TABLE]

where $\bar{\dot{n}}_{\rm SR}^{\rm iso}$ is the isotropic differential synchrotron photon production rate at each interaction point. The variables $V$, $t$, $\epsilon_{f}$ and $\Omega_{f}$ constitute the volume, time, outgoing photon energy, and outgoing photon solid angle, respectively,

For a power-law distribution of electrons with index $p$ between Lorentz factors $\gamma_{\rm min}$ and $\gamma_{\rm max}$ , denoted by the Heaviside step function $\Theta$, such that $\bar{n}_{e}(\gamma_{\rm e},\theta)=\Sigma_{\rm e}(\theta,\phi)/\Delta(\theta,\phi)(\gamma_{\rm e}/\gamma_{\rm min})^{-p}\allowbreak\times\Theta(\gamma_{\rm e}\;\gamma_{\rm min},\gamma_{\rm max})\Theta(\theta;0,\theta_{\rm max,X})$, the emissivity in the synchrotron power-law regime $\gamma_{\rm min}\ll\sqrt{\bar{\epsilon}_{f}B_{\rm cr}/\bar{B}_{s}}\ll\gamma_{\rm max}$ is given by Dermer & Menon (2009),

[TABLE]

where $\bar{B}_{s}(\theta,\phi)$ is the post-shock turbulent magnetic field magnitude in the comoving frame, and

[TABLE]

In the idealized MHD formulation where the bulk flow is force-free and there is no magnetic field perpendicular to the flow, we note that Lorentz transforming from the lab frame to the comoving blob frame preserves the parallel component $\langle B_{s}\rangle$.

To suitably develop light curves that would correspond to an observable, we form the differential synchrotron luminosity, $L_{\rm SR}(\Omega_{f},\epsilon_{f})\equiv d{\cal L}_{\rm SR}/(d\Omega_{f}d\epsilon_{f})$. This is the total source luminosity ${\cal L}_{\rm SR}$, or volume-integrated emissivity, binned in solid angle elements. We take advantage of the simple Lorentz transformation property of the spectral emissivity $\bar{j}_{\epsilon_{f}}$, namely that $j_{\epsilon_{f}}/{\epsilon_{f}^{2}}$ is invariant under boosts. Forming the differential luminosity in the observer frame, the result is an integral over the shock volume that is essentially a reconstruction of Eq. (10) of Dubus et al. (2015) in the optically thin limit:

[TABLE]

with the emissivity and photon energies being computed using the Doppler-shifted photon energy

[TABLE]

Observe that the factor $\delta_{\rm D}^{2}$ in Eq. (19), rather than $\delta_{\rm D}^{3}$ found routinely in relativistic jet contexts, arises because the integration over the emitting volume element is chosen to be in the observer’s frame rather than the comoving frame. This orbitally-modulated Doppler factor $\delta_{\rm D}$ is calculated for each point along the shock and dependent on the prescribed local bulk speed $\beta_{\Gamma}(\theta)$ and bulk Lorentz factor $\Gamma(\theta)$ along the shock, defined by

[TABLE]

where $\boldsymbol{\hat{u}}^{\prime}$ is the unit tangent vector along the polar direction of the shock in the inclined and orbital phase-rotated coordinate system (cf. Appendix A for definitions and conventions). Thus, we have

[TABLE]

In addition, $\boldsymbol{\hat{n}}_{v}$ is the unit vector in the direction of the observer, applicable to all orbital phases. Throughout, we employ the scaling Eq. (12)–(13), and associated $\beta_{\rm max}$ to prescribe the bulk flow speed along the shock.

The convolution of different Doppler factors realized in different portions of the integration volume near the shock is what defines the relative sharpness of peaks in the ensuing figures of orbital modulations of the flux. For the assumed symmetry, the integration Eq. (19) is two-dimensional in variables $\theta$ and $\phi$ for each point along the intrabinary shock. For a shock that is thin compared to the orbital length scale $a$, we may write the lab frame volume element as

[TABLE]

. The factor $\Delta(\theta,\phi)$ cancels with the factor in Eq. (17), and the net result is that the orbital-modulated differential luminosity or intensity is proportional to

[TABLE]

when $p$ is a constant along the region of interest. Crucially, the flux ratio of Eq. (25) generates energy-independent light curves when $p$ is spatially independent. In this regime, the ensuing large curves are largely regulated by geometric influences of the observer and ascribed velocity profile Eq. (12) . It is quickly discerned that the maximum amplitude of orbital modulation is attained when $\boldsymbol{\hat{n}}_{v}\boldsymbol{\cdot}\boldsymbol{\hat{u}}^{\prime}\approx 1/\beta_{\Gamma}$ realizing a flux enhancement of order $\lesssim 2\Gamma_{\rm max}^{2+(p-1)/2}$ modulo .

4.2 Application to PSR B1957+20, an SC-centered System

Using Eq. (4) we compute the volume-integrated emissivity ratio of superior-to-inferior conjunction in Figure 7, and take $\varsigma$ constant or $\varsigma(\theta)=\varsigma_{0}\left[1-(\theta/\theta_{\rm max,X})^{2}\right]$ at a fixed energy $\epsilon_{f}$ where the particle power-law is valid, with all numerical constant factors canceling. The ad hoc prescription of surface density profile $\varsigma(\theta)$ is known from hydrodynamic bow shocks, but likely does not trace the true time-averaged emitting particle distribution. The integration is taken to $\theta_{\rm max,X}=\pi/2$, that is, we only consider the head of the bow shock that participates for the hemisphere of the companion facing the pulsar. This portion of the shock geometry is expected to be largely axisymmetric and comprise the preponderance of accelerated charges and emission.

We employ the results of §2.2.1 to inform the choice of $R_{0}$ for a given inclination in the axisymmetric case. For lower maximum bulk speeds $\beta_{\rm max}=0.5$, Figure 7, although there is orbital modulation it is relatively modest and flat around SC. Little to no DP structure is exhibited except when the shock stagnation point is taken at the companion surface, assumed to be 90% of the Roche Lobe radius for B1957+20, where the influence of geometric occultation of the shock by the companion (solid curves) is larger. However, this DP morphology is rather flat and would require a large fraction of emission to be concentrated near the nose to produce the observed peak-dip ratio in Huang et al. (2012), . This small value of $R_{0}\approx R_{*}$ would also appear to be in tension with the radio eclipse estimates in §2.2.1 unless an optically-thin model for eclipses is operant. For other values of $R_{0}$, the influence of shadowing (solid versus dashed curves) is small especially for the more moderate inclination $i=65^{\circ}$.

Larger bulk speeds in the right panel in Figure 7 do produce characteristically DP synthetic light curves, which are in qualitative agreement with Huang et al. (2012) for the DP positions. As expected, modulation increases with inclination $i$ for either shock geometry. Given the large error bars on the data in Huang et al. (2012) and the early development stage of our model, we do not attempt any fitting. There is a complex interplay between the shock geometry, the system inclination, the value of $R_{0}$ (which principally moderates the influence of shadowing), the bulk Lorentz factor at points along the shock, and the surface density profile $\varsigma(\theta)$. Smaller inclinations reduce the peak separation for a given set of parameters, which can also be mimicked by altering the prescribed shock geometry for wider opening angles.

It is clear from this quantitative parameter exploration for B1957+20 that higher bulk Lorentz factors are preferred such that $\beta\approx 0.5-0.8$ is sampled by the observer line-of-sight cutting across the shock. This implies the shock components are not well-mixed unless the companion mass loss rate is lower than expected. Some concentration of emissivity near the stagnation point where Doppler-boosting is negligible, like the $\varsigma(\theta)$ prescription, is essential to mitigate the larger than observed amplitude and peak-dip ratio.

4.3 Application to PSR J1023+0038, an IC-centered System

In Figure 8 we compute the symmetric Doppler-boosted light curves with $p=1.2$ and $\theta_{\rm max,X}=\pi/2$, normalized to peak maximum, for a variety of inclinations and $\beta_{\rm max}=\{0.8,0.9,0.98\}$, the largest value corresponding to a maximum bulk Lorentz factor of about $5$. Since the shock surrounds the MSP, no shadowing by the companion is incorporated in the calculation since even a fully Roche lobe-filling companion for the given low inclinations $35^{\circ}-55^{\circ}$ only obscures a negligible fraction of the shock surface emissivity, and then so only near SC for emission that is Doppler de-boosted. As for B1957+20, we do not attempt any fits here, pending the development of a more complete and self-consistent model of particle transport, cooling and shock geometry, but rather explore shock geometry and physical parameters in a more generic way that will be useful to such a future study.

Finally, in Figure 9 we consider restricting $\theta_{\rm max,X}=\pi/6$, only considering emission from this small shock cap rather than the full head as in Figure 8. This scenario approximates one where emission from the shock is extraordinarily concentrated near the stagnation point. It is clear that such a restriction generally results in too-wide peak separation for any inclination or $R_{0}$ values, principally due to the local flatness of this region. Therefore, the emitting region in J1023+0038 and similar IC-centered systems must include synchrotron contributions from larger $\theta_{\rm max,X}$ values where the shock significantly curves or bows, consistent with large radio eclipse fractions enshrouding the MSP.

5 Conclusion

In this paper we have focused on the radio and X-ray phenomenology of MSP binaries constructing geometric models for X-ray double-peaked light curves and radio eclipses. We have shown that in the thin-shell model, one may constrain the shock parameters as a function of binary inclination using radio eclipses and X-ray light curves.

We have argued that the shock in BWs and RBs have two components that are not well-mixed unless the companion mass loss rate is low, consistent with the reasoning that relatively low baryon loading is necessary for the mildly relativistic bulk Lorentz factors attained in the shocked pulsar wind component that generates the DP Doppler-boosted light curves. The inferior conjunction phase-centered double-peaked light curves imply a relatively stable shock at orbital length scales in our model. Accordingly, some mechanism is necessary to stabilize the shock from gravitational influences for these cases to prevent the systems from readily transitioning to an LMXB state.

The geometric explorations of the shock presented here will be used in a future paper for modeling high-energy emission and transport processes. The quantitative agreement of the synthetic light curves with observations strongly encourages further development of the model to include considerations of diffusive and convective transport in the shock environs. There are implications for orbitally-modulated inverse Compton emission in the Fermi LAT and TeV bands, depending on the shock geometry and target photons from the shock and companion. Model fitting explorations in the near future of DP X-ray light curves may be able constrain the shock physics given orbital parameters, or constrain orbital parameters such as binary inclination (and consequently MSP mass) using a model for the shock geometry and the particle distribution in different regions. Energy-dependent light curves by NICER or NuSTAR may also probe the underlying relativistic shock acceleration’s spatial dependence in the future. Our study lastly motivates population studies of radio eclipse phenomenology in the MeerKAT/SKA era, where the number of MSPs discovered will increase by severalfold (Keane et al., 2015).

We thank the anonymous referee for aiding to improve the flow and structure of the manuscript. Z.W. acknowledges helpful discussions with Cees Bassa, Julia Deneva, Guillaume Dubus, Jason Hessels, Christopher Johns-Krull, , Edison Liang, Alessandro Papitto, , Scott Ransom, Mallory Roberts, Bronek Rudak, Ben Stappers, and Kent Wood. C.V. & Z.W. are supported by the South African National Research Foundation (NRF). The work of M.B. is supported by the South African Research Chairs Initiative of the Department of Science and Technology and the NRF. Any opinion, finding and conclusion or recommendation expressed in this material is that of the authors and the NRF does not accept any liability in this regard. A.K.H. and M.G.B. acknowledge support from the NASA Astrophysics Theory Program. A.K.H., Z.W., and C.V. also acknowledge support from the Fermi Guest Investigator Cycle 8 Grant.

Appendix A Eclipses by Surfaces in Binary Systems

A.1 General Formalism

Although the focus of this text is on eclipses by optically-thick intrabinary shocks, the analytical formalism here can also be applied to arbitrary azimuthally symmetric surfaces that occult a point source. The case of eclipses by quasi-static Roche lobes has been treated analytically in Chanan et al. (1976) and Kopal (1959), however the method presented here is more broadly applicable.

Without loss of generality, we prescribe a right-handed orthonormal cartesian coordinate system, in flat spacetime, such that the barycenter is the origin and the orbital angular momentum vector of the binary is in the $\boldsymbol{\hat{z}}$ direction with binary angular frequency $\Omega_{\rm b}$. The observer angle is defined such that $\boldsymbol{\hat{z}\cdot\hat{n}_{v}}=\cos i$ with $\boldsymbol{\hat{n}_{v}}=\cos i\boldsymbol{\hat{z}}+\sin i\boldsymbol{\hat{x}}$ the observer line-of-sight unit vector. To obtain the projection of the binary system into the plane of the sky perpendicular to $\boldsymbol{\hat{n}_{v}}$, we construct a new primed rotated coordinate system such that $\hat{\bf{x^{\prime}}}$ is parallel to $\hat{\bf{n}}_{v}$ at arbitrary orbital phase. If the cartesian basis $\boldsymbol{\hat{r}}\rightarrow(\hat{\bf{x}},\hat{\bf{y}},\hat{\bf{z}})^{\intercal}$ defines the vector space at phase $\Omega_{\rm b}t=0$ and inclination $i=\pi/2$, then the primed orthonormal coordinate basis for arbitrary orbital phase and inclination is constructed by two successive rotations. Here the inclination is formulated such that $i=0$ and $i=\pi/2$ are face-on and edge-on views, respectively, and with the phase convention $\Omega_{b}t=0$ prescribing the superior conjunction where the companion is between the MSP and observer in Eq. (A3) of the orbital equations of motion. The rotations are given by

[TABLE]

where

[TABLE]

define the primed coordinate basis about the barycenter such that $\hat{\bf{x^{\prime}}}=\hat{\bf{n}}_{v}$, with the span of $\hat{\bf{y^{\prime}}}$ and $\hat{\bf{z^{\prime}}}$ characterizing the plane of the sky.

The stars are treated as following unperturbed Keplerian trajectories, although a post-Newtonian generalization is straightforward, with phase zero constructed to rest on the $\hat{\bf{x}}$-axis with the following equations of motion,

[TABLE]

The scalars $r_{i}(t)$ for non-zero eccentricity depend on the orbital phase,

[TABLE]

where $r_{\rm c}+r_{\rm NS}\equiv 1$ (normalized to $a$) are the respective semi-major axes for each star, with $r_{\rm c}=q/(q+1)$ and $r_{\rm NS}=1/(q+1)$ where $q=M_{\rm MSP}/M_{\rm c}$ is the mass ratio. For nonzero eccentricity, which is beyond the scope of this paper, an additional rotation of coordinates by the argument of periastron must be incorporated in Eq. (A1).

Rotating coordinates to the primed coordinate basis by Eq. (A1), the position vectors of the primary and secondary are given by

[TABLE]

with a parallel/orthographic projection on the ${\rm Span}\{\hat{\bf{y^{\prime}}},\hat{\bf{z^{\prime}}}\}$ plane comprising the observer view.

For orbital angular speeds $a\,\Omega_{\rm b}\ll c$, the coordinates of a vector $\boldsymbol{\ell}$ defining a ray leaving the primary towards the observer traversing a distance $l$ (in units of $a$) can be expressed as

[TABLE]

with the distance between an interaction point at the ray and an arbitary location in the system $\boldsymbol{r}_{\rm m}^{\prime}$ (e.g., the secondary) given by $\left|\boldsymbol{\ell}^{\prime}-\boldsymbol{r}_{\rm m}^{\prime}\right|$. The optical depth for a specified absorption coefficient $\alpha$ is computable in the usual way through the observer line-of-sight integral at each orbital phase, $\tau=\int\alpha dl$. Such a procedure for pulsar eclipses by specified volumes has been performed in the context of MSP binaries (Rasio et al., 1989) and other pulsar systems (e.g., Lyutikov & Thompson, 2005) but relies on a physical model underpinning the absorption coefficient, and is computationally inefficient for arbitrary volumes when the optical depth is large for the region of interest, and when the transition region of low-to-high optical depth is sharp. In the case of optically thick surfaces and volumes, we employ a geometric occultation formalism described below.

For a surface with azimuthal symmetry with radial function $R(\theta)$, every locus of points at fixed $\theta$ is a circle. When projected onto the ${\rm Span}\{\hat{\bf{y}}^{\prime},\hat{\bf{z}}^{\prime}\}$ plane of the sky, these circles are transformed to ellipses by the elementary rotations of orbital phase and inclination angle. The definition of eclipses by a bow shock is thus transformed to the problem of determining if the projected location of the MSP is interior to any projected bow shock ellipse, for all prescribed $\theta$. We define $\theta_{\rm max,R}$ as the maximum polar angle where the shock surface transitions from optically-thick to thin at a given frequency. The maximum length along the axis of symmetry of the shock tail past the companion position $L_{z}(\theta_{\rm max,R})=-R(\theta_{\rm max,R})\cos\theta_{\rm max,R}$ should be smaller than the typical orbital length scale $a$ for flow velocities that are similar in scale to the orbital speed.

The locus of points for a particular projected ellipse, given the conventions developed above, can be found to be

[TABLE]

with

[TABLE]

where $\phi\in(0,2\pi)$ is the azimuthal angle parameterizing the azimuthally-symmetric surface. The function $S(t,\theta)$ defines whether the shock surface surrounds the primary or secondary with the orbital phase convention defined above: for surrounding (bowing around) the secondary it takes the form $S_{\rm c}=r_{\rm c}(t)-R(\theta)\cos\theta$ but $S_{\rm NS}=-r_{\rm NS}(t)+R(\theta)\cos\theta$ for where the shock enshrouds the pulsar.

With these definitions, it is evident that the center of the ellipse is then $\{\alpha_{y},\alpha_{z}\}$. In general, the ellipses’ axes of symmetry are rotated by an angle with respect to the $\hat{\bf{y}}^{\prime}$ or $\hat{\bf{z}}^{\prime}$ axes. The semi-major $E_{a}$ and semi-minor $E_{b}$ axes lengths are easily shown to be

[TABLE]

for a given inclination and orbital phase and, for a given $\theta$. The projected eclipsed area is $\propto R(\theta)^{2}\sin^{2}\theta\sin i\cos(\Omega_{\rm b}t)$. Using elementary methods, the angle between an ellipse’s semi-major axis and the $\hat{\bf{y}}^{\prime}$ direction is found to be

[TABLE]

for phases $-\pi/2\leq\Omega_{\rm b}t\leq\pi/2$.

These relations then allow us to analytically test whether the MSP’s projected coordinate is eclipsed by an ellipse. Rotating and translating coordinates, we construct new axes aligned with the ellipse with origin at the ellipse center. The coordinates of the MSP are then given by

[TABLE]

which can be used to test if the projected position of the MSP is eclipsed. We define a unit Heaviside test function $\Theta(\Gamma)$ with argument

[TABLE]

is less than zero for eclipses. Hence the position $G_{\rm x,y}$ and parameter $\Gamma$ can be associated with a physical attenuation (and $\tau$) at any given $\theta$, i.e., $e^{-\tau}\sim f(G_{\rm x},G_{\rm y})$ for some function $f$; we choose $f$ to be a unit step function in our optically thick formalism. Then, $\Theta(\Gamma)$ is unity when the pulsar is not eclipsed. The eclipsing by the entire ensemble of ellipses $j_{\rm max}$ that encompass the employed surface is a product of $\Theta$ for every $\theta$ up to $\theta_{\rm max,R}$, suitably discretized to sample the complete surface,

[TABLE]

For the case where the surface surrounds the emission point source and is optically-thick, a lower limit of the shrouding fraction can be found by considering the condition Eq. (A19) evaluated only at $\theta_{\rm max,R}$. This yields a lower limit because a general surface may have pockets of plasma that may occlude the emission source in a nontrivial manner. For cases where the bow shock is not swept back and the curvature is positive everywhere, the preceding construction is unnecessary, since the shrouding fraction is analytically expressible by Eq. (11).

It is numerically expedient to adopt a continuous approximation of the Heaviside function, $\Theta[\Gamma(\theta_{j})]\approx 1/(1+e^{-{\cal T}\Gamma(\theta_{j})})$, where ${\cal T}$ is an associated pseudo optical depth that parametrizes the geometric transition thickness near the occluding surface $R(\theta)$ and approximates a sharp transition for ${\cal T}\gg 1$ for a fixed observational frequency. This transforms the product to a sum,

[TABLE]

Because of the symmetry in the problem, the full orbit and other branches for the solution of $\cos\vartheta_{E}$ need not be considered, and the calculation can be restricted to one quarter of the total orbit. The total eclipse fraction $f_{E}$ by the surface during a total duration $2\pi$ of an orbit is then given by the integral,

[TABLE]

for a given set of Keplerian orbital parameters and radial function $R(\theta)$.

A.2 “Cometary” Tail Sweepback Due to Orbital Motion

If the locus of points of the shock flow is azimuthally symmetric relative to the companion, i.e., axisymmetric at each instantaneous orbital phase then it is straightforward to include the sweep-back and Coriolis forces due to orbital motion by retarding the surface at each $\theta$ by a specified time delay. We neglect Coriolis effects perpendicular to the outward direction, equivalent to the assumption that $v_{\perp}\ll v_{\rm inj}$ where $v_{\perp}$ is the bulk flow velocity perpendicular to the line connecting the two stars; this approximation, which preserves axisymmetry of the shock, is generally valid for locales distant from the stagnation point, the principal focus of this work associated with radio eclipses. Connecting orbital phases to the surface parameter $\theta$ then suffices to produce a good approximation of such swept-back shock tails, and the eclipse fraction formalism of the previous section can be expediently utilized. This ballistic assumption is valid if the shock bulk velocity is highly supersonic, as in MSP binaries, especially if the winding number of the swept-back tail that contributes to the eclipses is zero.

The model we adopt here simply injects the shock flow at a finite velocity $v_{\rm inj}$ near the companion at $\theta=\pi/2$, assumed much greater than the escape velocity $v_{\rm inj}\gg v_{\rm esc}\sim 10^{7}$ cm s*-1* which, coincidently, is also the same order of magnitude as the orbital velocity $v_{\rm orb}$, so that gravitational effects can be neglected for BWs. More complex self-consistent prescriptions for the bulk flow can also adopted in this framework, but may introduce more parameters. The time delay for nonrelativistic velocities and finite constant acceleration $a_{\rm wind}$ (in the lab nonrotating frame) due to radiation pressure from the pulsar wind is

[TABLE]

which reduces to $aL_{z}/v_{\rm inj}$ for small flow accelerations. For zero acceleration, the approximate tail angle with respect to the line connecting the two stars is then $\arctan(v_{\rm orb}/v_{\rm inj})$, good when $v_{\rm inj}\gg v_{\rm orb}$. For large $a_{\rm wind}$ or $v_{\rm inj}$, the time delay approaches zero and the symmetric shock solution is recovered. Here $L_{z}(\theta)=-R(\theta)\cos\theta$ connects the time delay to points along the shock surface. Thus, generically, increasing asymmetry increases the total eclipse fraction at a given $\theta_{\rm max,R}$.

The computation of eclipse fraction is straightforward with the simple replacement $\Omega_{\rm b}t\rightarrow\Omega_{\rm b}(t^{\prime}-\delta t)$ in Eq. (A18) at each $\theta>\pi/2$ up to a maximum $L_{z}$ or $\theta_{\rm max,R}$ for all variables associated with the shock surface. The ingress and egress of eclipses are defined to occur for phases $\Omega_{\rm b}t^{\prime}\in[-\pi/2,0)$ and $\Omega_{\rm b}t^{\prime}\in(0,\pi/2]$ respectively. The total eclipse fraction is then $f_{E}=f_{E}({\rm Ingress})+f_{E}({\rm Egress})$, and the asymmetry can be characterized by the difference $\Delta f_{E}\equiv f_{E}({\rm Egress})-f_{E}({\rm Ingress})$ and the ratio of egress to ingress eclipse durations.