Constraining Relativistic Bow Shock Properties in Rotation-Powered Millisecond Pulsar Binaries
Zorawar Wadiasingh, Alice K. Harding, Christo Venter, Markus Böttcher, Matthew G. Baring
DDraft version July 18, 2018
Typeset using L A TEX default style in AASTeX61
CONSTRAINING RELATIVISTIC BOW SHOCK PROPERTIES IN ROTATION-POWERED MILLISECONDPULSAR BINARIES
Zorawar Wadiasingh, Alice K. Harding, Christo Venter, Markus B¨ottcher, and Matthew G. Baring Centre for Space Research, North-West University, Potchefstroom, South Africa Astrophysics Science Division, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA Department of Physics and Astronomy, Rice University, Houston, TX 77251, USA
ABSTRACTMultiwavelength 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 interpretthe radio to X-ray phenomenology in a consistent framework. We advocate the existence of two distinct modesdiffering in their intrabinary shock orientation, distinguished by the phase-centering of the double-peaked X-ray orbitalmodulation 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 . We develop synthetic X-raysynchrotron orbital light curves and explore the model parameter space allowed by radio eclipse constraints appliedon archetypal systems B1957+20 and J1023+0038. For B1957+20, from radio eclipses the stand-off is R ∼ .
15– 0 . R (cid:46) . . (cid:46) R (cid:46) . Keywords: radiation mechanisms: non-thermal — pulsars: individual (J1023+0038, B1957+20) —binaries: eclipsing — X-rays: binaries
Corresponding author: Zorawar [email protected] a r X i v : . [ a s t r o - ph . H E ] M a r Wadiasingh et al. INTRODUCTIONThe old population of rapidly-spinning neutron stars, generally known as the millisecond pulsars (MSPs), are fre-quently found as binaries. In the standard “recycling” evolutionary scenario, MSPs attain their short rotation periodsthrough angular momentum transfer by accretion from a main-sequence companion in a low-mass X-ray binary phase(Alpar et al. 1982). Depending on the initial conditions, such evolution can yield an MSP binary with low-masscompanion ( (cid:28) . M (cid:12) ) in a circular orbit with a short orbital period < γ -rayMSP binaries in tight circular orbits with low-mass companions are useful astrophysical laboratories for the physicsof pulsar winds and relativistic shock acceleration. They are relevant to not only striped pulsar winds but also tothe 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 toobserver-dependence of sampling, via orbital modulations, of the emission region in a viewing geometry constrained byradio and optical determinations of the binary mass functions. Prior to the launch of the Fermi
Large Area Telescope(LAT; Atwood et al. 2009) γ -ray observatory, only three such “black widow” low-mass radio MSP binaries were knownin the Galactic Field. The first of these is the original “black widow” B1957+20 (Fruchter et al. 1988), a 1 . M c (cid:38) . M (cid:12) with a binary period of 9 .
17 hours. It is now well-establishedthat old “recycled” MSPs emit γ -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 starand its wind. Accelerated leptons from MSP binaries may also significantly contribute to the anomalous rise in theGalactic energetic positron fraction observed in low-Earth orbit (Venter et al. 2015) by the Alpha Magnetic Spectrom-eter (AMS-02), PAMELA , and
Fermi
LAT instruments. Follow-up observations, predominantly in the radio band, of
Fermi
LAT unidentified sources have expanded the binary population to over 30 in the Field, bringing their totalnumber to over 70 known when including those residing in globular clusters (Manchester et al. 2005) with the caveatthat the companion mass is unclear in many of these systems.Precision radio timing of the binary MSPs accounting for orbital Doppler wobbles yields the pulsar binary massfunction andan estimate of the minimum mass of the companion and semi-major axis of the orbit, typically ∼ cm,based on inclination sin i ≤ . M (cid:12) .Empirically, the known population of rotation-powered MSP low-mass short-period binaries are loosely segregatedbased on the minimum companion mass M c (Roberts 2011): black widows (BWs) with minimum companion masses M c (cid:46) . M (cid:12) that may be degenerate, and the rarer redbacks (RBs) with non-degenerate companions M c (cid:38) . M (cid:12) .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 pulsarwind. Recently, some RBs have been observed to transition between a rotation-powered pulsar state and a low-massX-ray binary accretion state (Archibald et al. 2009; Papitto et al. 2013), confirming the association between thesesource classes. The behavior of these transitional objects is complex, exhibiting poorly understood transitions inX-ray luminosity and accretion states (Linares 2014; Bogdanov et al. 2015). The focus of this paper is modeling theconceptually clearer rotation-powered state of BWs and RBs where the total energy budget is constrained by thepulsar rotational energy loss rate ˙ E SD .In the standard narrative envisioned for BWs and RBs, an intense ˙ E SD ∼ − erg s − MSP pair plasma“striped” wind reprocesses in an intrabinary shock that irradiates the tidally-locked companion, preferentially heatingthe facing side. The companion exhibits orbitally-modulated optical light curves interpreted as a convolution ofellipsoidal variations due to the companion being nearly Roche lobe-filling and anisotropic photospheric emission dueto pulsar irradiation. Besides heating by irradiation, particle acceleration beyond TeV energies can take place in therelativistic magnetized intrabinary shock (Harding & Gaisser 1990; Arons & Tavani 1993), whose pressure supportderives from either the companion wind or a magnetosphere. The companion wind matter also generates radio eclipsesof the MSP at orbital phases where the companion is between the pulsar and observer. In this picture, eclipsesin the radio pulsations of the MSP in some systems can assist in constraining the shock and wind geometry as wedemonstrate in this paper. A photometric light curve model of the orbital optical variations of the companion can beused to constrain the system inclination and irradiation efficiency (e.g., Breton et al. 2013). If radial velocities of the https://confluence.slac.stanford.edu/display/GLAMCOG/Public+List+of+LAT-Detected+Gamma-Ray+Pulsars hocks in MSP Binaries M MSP well-above the canonical 1 . M (cid:12) (e.g., vanKerkwijk et al. 2011). Such massive stars constrain theories of the nuclear equation of state, marking these systemsas attractive targets for the Neutron star Interior Composition Explorer (NICER, Arzoumanian et al. 2014) due forlaunch in 2017 with an energy range of 0 . XMM-Newton .The scrutiny of BWs and RBs can help uncover the largely unknown physics of pulsar winds in MSPs, and aid inunderstanding where the transition occurs from a magnetic flux-dominated to a particle-dominated flow. Unlike theCrab Nebula whose termination shock or inner knot is ∼ − cm away from the pulsar, the companions in BWsprovide a fixed target at a distance only ∼ cm from the pulsar with much higher magnetic fields realized thanin PWNe shocks. Indeed, the “clean” nature of the circular orbit, tidally-locked companion, steady well-constrainedpulsar spin-down energy budget, and multiwavelength observations establish these systems as useful probes for studyingthe physics of pulsar winds and shock acceleration. 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 processfor leptons if the pulsar wind magnetization σ , the ratio of magnetic to pair plasma particle kinetic energy density, islarger than unity – however too large a σ may preclude the existence of the observed shock. Conversely, if σ is small, amore conventional diffusive shock acceleration (DSA) may be the energization mechanism but is likely less efficient dueto 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 mightbe 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 ratherthan the pulsar since the pulsar wind ram pressure dominates that of the companion wind. Although many physicalprocesses 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 and elaborated in § R LC = cP MSP / (2 π ), where P 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, with gravitationalinfluences of the MSP significantly affecting the companion wind. The shock geometry may significantly impact modelsof 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 analysesof 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 andRBs (see § Fermi
LAT and the planned ˇCerenkov Telescope Array (CTA). In § §
3. In § § CONSTRAINING THE INTRABINARY SHOCK IN BLACK WIDOW AND REDBACK SYSTEMS2.1.
Overview and Idealizations
As is generically the case for collisionless astrophysical shocks, there are two components separated by the contactdiscontinuity for the intrabinary shock in MSP binaries like B1957+20 and J1023+0038: a relativistic pair shockedpulsar wind and an ionized shocked companion component (Phinney et al. 1988; Eichler & Usov 1993). Such a gross
Wadiasingh et al.
Figure 1.
Schematic cross-sectional diagram of a “spider” MSP binary system, scale exaggerated for clarity, illustratinggeometry and defining some variables and parameters used in this paper. The gray dotted curves depict differing thin-shellshock surface realizations for colliding momentum-dominated isotropic winds. In RBs and transitional objects in the rotation-powered state, the orientation of the shock is reversed such that the shock bows around the pulsar somewhat outside its lightcylinder, with the definition of R taken from the MSP in that case. Definitions of symbols are discussed in § structure must exist and is borne out in hydrodynamic (Bucciantini 2002; van der Swaluw et al. 2003; Bosch-Ramonet al. 2012; Dubus et al. 2015) and relativistic magnetohydrodynamic (RMHD, e.g., Bucciantini et al. 2005) simulationsof pulsar wind shocks in various contexts. Relativistic plasma/magnetic turbulence and electron acceleration, likelymediated by magnetic reconnection, DSA, or energization mediated by shear flows (Liang et al. 2013) if mixing betweencomponents 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 parallelshocks and may lead to spatial dependence in the acceleration in the bowed “head” of the intrabinary shock. Due tothe disparate length scales of the shock and gyroscale acceleration, such particle acceleration by current kinetic-scalesimulations (e.g., particle-in-cell codes) cannot be computed in a self-consistent manner over the large length scalesof the shock. Developing an expedient formulation to empirically diagnose the spatial character of such accelerationfrom high-energy spectroscopically phase-resolved light curves is a planned goal for future work.A simplified structure of the pulsar binary and intrabinary shock, central to this paper, is depicted in Figure 1. Thecircular binary components orbit with radii r c or r NS around a common center-of-mass, with separation a = r c + r NS for a mass ratio q = M MSP /M c (cid:29) R ∗ that is ≤ L Lagrange point distance. This spherical approximation for the companion, adopted for expediency, is employed forshadowing and eclipsing calculations in § § ˆΩ b is inclined at angle i withrespect to the observer line-of-sight ˆ n v ; if the pulsar spin axis is aligned with ˆΩ b as one may expect from the recyclingevolution, then ζ = i . The MSP’s pulsed radio emission, which originates within the light cylinder, is assumed to bepoint-like. This is a good approximation since the pulsar magnetosphere is small compared to the orbital separation a , i.e. R LC /a ∼ − for a ∼ cm and a typical MSP spin period of 2 ms. This small distance scale relative to hocks in MSP Binaries a is also approximately the MSP striped relativistic MHD wind wavelength length scale (all distances hereafter arespecified in units where a = 1 unless otherwise noted). Contours representing the intrabinary shock surface in thethin-shell approximation for two colliding isotropic momentum-dominated winds (Canto et al. 1996) are shown, withthe purple curve highlighting a particular case at the stagnation point R . The polar angle θ max , R defines the regionthat is optically thick at a particular radio frequency for radio eclipses of the MSP in § § §
4. The relevance of red labels referring to fast/slow cooling willbe apparent in due course in § i may be constrained using severalmethods: radio MSP and optical companion mass functions, orbitally-modulated γ -ray emission and eclipses, radioeclipses, and orbitally-modulated intrabinary X-ray shock emission. The latter two constitute the purview of thispaper where i is a critical model parameter. In addition, if ζ ≈ i , pulsed radio and γ -ray light curve models (e.g.,Guillemot & Tauris 2014; Johnson et al. 2014) and pulsed thermal X-rays from polar-cap hot spots of the MSP informon the orbital geometry.Phase zero of the orbit in most observational contexts for MSP binaries is defined where the pulsar passes theascending node, with orbital phases 0 .
25 and 0 .
75 the superior conjunction (SC) and inferior conjunction (IC) of theMSP, respectively — in this article, phase zero is defined as SC where the MSP is behind the companion for theobserver as this is a natural choice for radio eclipses. For some MSP binaries where the pulsar is totally shroudedin the radio and no radio ephemeris is available, the IC phase is associated with the global optical maximum of theirradiated 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 syn-chrotron cooling of relativistic electrons and positrons at an intrabinary shock in a turbulent and relatively high (cid:46)
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 emis-sion than for dipping/eclipsing LMXBs (e.g., Parmar et al. 1986) where accretion power may dominate. The emissiontypically has a strong non-thermal power-law component, with relatively flat photon indices Γ X ∼ − . p ≈ X − (cid:46) ∼ ∼ < ∼ Wadiasingh et al.
Table 1.
Rotation-Powered BWs and RBs with DP X-Ray Light CurveMorphologyName 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]
Note —Current list of MSP binaries in the rotation-powered state forwhich DP X-ray emission attributed to Doppler-boosting has been ob-served. IC and SC denote inferior and superior conjunction of the pulsar,respectively.
References — [1] Stappers et al. (2003); Huang et al. (2012) [2] Bogdanovet al. (2005) [3] Archibald et al. (2010); Tam et al. (2010); Bogdanovet al. (2011, 2014b); Tendulkar et al. (2014) [4] Bogdanov et al. (2014b);de Martino et al. (2015) [5] Bogdanov et al. (2014a); Hui et al. (2014) [6]Romani (2015); Salvetti et al. (2015) [7] Roberts et al. (2015); Hui et al.(2015) [8] Gentile et al. (2014); Romani et al. (2015) [9] Romani & Shaw(2011); Yatsu et al. (2015) the DP morphology, typically ranges 65 −
85% of the global maximum and may be at some small phase-offset fromIC or SC. Although there is some statistical uncertainty, peaks are generally not identical with the leading peak oftenmore prominent than the trailing peak. Peak-to-peak separation is confined to a rather narrow range of 0 . − .
35 innormalized phase for the sources in Table 1 while peak full-widths are generally around 0 .
1. Interestingly, J2129–0429which exhibits one of the most well-defined DP morphologies is also close to edge-on with i ≈ ◦ (Bellm et al. 2016).Simple occultation by the companion of the emission region 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) andTendulkar et al. (2014), since this is 0 . ◦ to 53 ◦ 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 radioshrouding of the MSP suggests the plasma is not well-confined near the companion. However it is difficult to envisiona plausible and relatively stable hydrodynamic scenario where a shock exists near the companion L point but otherplasma is shrouding the pulsar (cid:38)
50% of the orbit but generally not at pulsar IC for such low (cid:46) ◦ inclinations.Moreover, for an X-ray emission region close to the companion, occlusion also innately leads to a DP structure thathas a peak separation of ∼ . the phase centering of the DP structure is a key discriminant of the shock orientation and system state . Inaddition, the light curve structure serves as a probe of shock geometry, particle acceleration, and shock mixing. hocks in MSP Binaries + e − wind and the shock-heated ionized companion matter, that is, the baryon loading of the flow.For a striped wind of magnetization σ 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 fieldorientation, the shock geometry is quasi-perpendicular at the nose with the highest compression ratio, transformingsmoothly to quasi-parallel at the flanges with a lower compression ratio. This spatial dependence of the compressionratio, relativistic shock obliquity, along with higher particle resident time near the stagnation point (the fast coolinglocale in Figure 1), should inherently influence the local particle acceleration, cooling and emergent radiation dependingon what shock locales the observer line-of-sight samples as a function of orbital phase. However, a detailed explorationin a self-consistent geometry with a transport model for leptons along the shock is deferred to a future paper. For ourpresent study, we focus on the gross DP structure of the light curves in different geometries that can easily be adaptedfor different sources and energies.It can be shown that the equatorial upstream wind magnetic field magnitude B w , dominated by the toroidal com-ponent at large cylindrical radii r s (cid:29) R LC from the pulsar, is B w ≈ (cid:32) E SD c (cid:33) / r s = 22 (cid:32) ˙ E SD erg s − (cid:33) / (cid:18) cm r s (cid:19) G . (1)This relatively large magnetic field advocates synchrotron cooling as a significant energy loss mechanism for electrons.A rudimentary estimate for the pulsar contribution to the electron/positron number density near the shock may befound by assuming isotropic particle outflow from the MSP at a multiplicity M ± of the Goldreich-Julian rate ˙ N GJ (Goldreich & Julian 1969) from the pulsar polar caps,˙ N GJ ≈ cA cap | ρ GJ | e ≈ √ ce ˙ E / s − , (2)where | ρ GJ | = | Ω · B | / (2 πc ) ∼ B/ ( cP MSP ) the Goldreich-Julian charge density and A cap ≈ πR (1 − (cid:112) − R MSP /R LC ) is the approximate pulsar polar cap area for an aligned rotator. Then for a secondary pairmultiplicity M ± , the pulsar contribution to the number density at distance 10 cm is n e , MSP = M ± ˙ N GJ (4 πcr s ) ≈ × − M ± (cid:32) ˙ E SD erg s − (cid:33) / (cid:18) cm r s (cid:19) cm − . (3)For MSPs, the secondary multiplicity from pair cascade codes is typically M ± ∼ – 10 of the primary polar capoutflow rate (Harding & Muslimov 2011; Timokhin & Harding 2015; Venter et al. 2015) while constraints from youngPWNe studies (Sefako & de Jager 2003) or the Double Pulsar (Breton et al. 2012) suggest M ± ∼ − . Thus forBWs, the typical pulsar contribution probably does not exceed ∼ cm − unless the pair wind is highly anisotropicin the plane of the orbit. For IC-centered spiders and transitional systems where the shock may be much closer tothe MSP, the pulsar pair density can be profoundly larger by a factor up to ( a/R LC ) (cid:46) and may be a significantinfluence for the radio eclipses and radiation physics.For a well-defined MHD shock to develop, the magnetization must attain σ (cid:28) σ = B / (4 π n e , MSP (cid:104) γ w (cid:105) m e c ) (cid:46) (cid:104) γ w (cid:105) for an isotropic pair wind, (cid:104) γ w (cid:105) (cid:38) (cid:32) E SD c (cid:33) / e M ± m e c ≈ × M ± (cid:32) ˙ E SD erg s − (cid:33) / , (4)where M ± (cid:29) σ (cid:28)
1, themagnetic field in the shocked pulsar wind field B s then scales as B s ∼ √ σB w in the ultrarelativistic perpendicularshock limit (Kennel & Coroniti 1984). However, the magnetic dissipation processes upstream may convert or destroythe striped wind morphology such that the shock may be quasi-parallel in the proper frame. A containment argument, Wadiasingh et al. based on the observed X-ray power law provides a rudimentary lower bound on B s – the Larmor radius r L of electronsin the shock must be smaller than about 1% of the orbital length scale r L (cid:46) . a ∼ cm. Then, assuming emissionat the critical synchrotron dimensionless energy (cid:15) c = 3 B s / (2 B cr ) γ with B cr ≈ . × G and electron Lorentzfactor γ e , for an observed power law extending to energy (cid:15) X, max in units of m e c , B s (cid:38) B s , min ≈ . (cid:15) / X, max (cid:18) cm r L (cid:19) / G , (5)where we have neglected factors of roughly unity associated with Doppler shift of energies corresponding to mildlyrelativistic bulk speeds along the shock. Therefore power laws extending up to (cid:15) X, max ≈ . ≈
80 keV / ( m e c )observed by NuSTAR for J1023+0038 advance 2 (cid:46) B s (cid:46) B w ∼
200 G in the relativistic magnetized shock if r s ∼ cm, which implies radiating electron Lorentz factors of order 10 –10 , i.e. well-above a thermal population. A moreloose assumption of r L ∼ a still results in B s (cid:38) − G, still considerably higher than those in PWNe. Thereforethis synchrotron component extends into the UV/optical/IR and lower energies, but such a power-law extrapolationyields expected fluxes well-below the sensitivity of any facility. For other spiders where observations at energies abovethe classical soft X-ray band are not available, the field magnitude is still greater than about one Gauss, orders ofmagnitude larger than those in plerions. We consider implications of these bounds on the shock in § Radio Phenomenology
Orbital eclipses of the MSP’s radio pulsations are a common feature in many BWs and RBs in the rotation-poweredstate. Observed orbital eclipse fractions f E are ordinarily f E ∼ −
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 ∼
60% of an orbit at 150 MHz (Archibald et al. 2009, 2013). Some BWs also have extensive eclipses. There appearsto be a dichotomy in the relative stability of eclipses – for some BWs like B1957+20 eclipses near SC are generallystable orbit-to-orbit, while sporadic mini-eclipses are seen in some other systems particularly those systems with largereclipse 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 symmetricand antisymmetric parts about SC is attainable as a function of observer frequency ν . Frequency dependence of theeclipse 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). Atthe highest radio frequencies ν , 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 aresignificantly larger than R ∗ for a fully Roche lobe-filled companion, even for sin i ≈
1. No eclipses by the companion areexpected if i < ◦ − arcsin( R ∗ /a ), but many systems with eclipses have well-constrained inclinations and companionsizes which violate this inequality. Therefore, eclipses must be predicated on plasma within the system and/or asecondary magnetosphere. Eclipses typically exhibit large plasma dispersion measures before the coherence in thetiming solution of pulsations is lost, likely due to absorption rather than scattering (Roy et al. 2015); continuumeclipses 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 ∝ ν − . .There are a panoply of potential eclipse mechanisms (cf. Michel 1989; Eichler 1991; Gedalin & Eichler 1993; Thompsonet al. 1994) depending on physical parameters realized in the intervening plasma. Cyclotron absorption has beenposited in B1957+20 (Khechinashvili et al. 2000) but relatively little Faraday rotation is seen, consistent with a 1 −
10G mean magnetic field magnitude in the eclipsing medium (Fruchter et al. 1990), not inconsistent with Eq. (5) sincethe eclipsing medium consists of the ionized companion wind as well. Moreover, it is now known that the companion inB1957+20 is likely non-degenerate (Reynolds et al. 2007). Excess delays, consistent with plasma dispersion, generallyshow that the average free electron column density rises sharply from (cid:104) n e (cid:105) d ∼ cm − to 10 cm − at phasesdeep into the eclipse (Ryba & Taylor 1991; Stappers et al. 2001) for BWs, for d ∼ a the line-of-sight column depth,but it is anticipated that there is also clumping near the shock contact discontinuity. This (cid:104) n e (cid:105) is much higher thanimplied by Eq. (3), 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 geometricvolume 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 surfacegermane to the X-ray light-curve phasing in Table 1. For the SC-centered DP phase centering where the shock is hocks in MSP Binaries 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 erraticeclipses are expected where the companion wind can enshroud the pulsar, and is necessarily turbulent for the obligatoryangular momentum loss. The radio optical depth, as well as the shock orientation depend on the companion windmass 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 windsource since the circularization radius R circ must be larger than the shock radius R (measured from the MSP), or thesystem will be predisposed to a disk-state (Frank et al. 2002). Moreover, for the radio pulsar state, R must exceedthe light cylinder scale, that is, R > Max( R circ , R LC ). This then favors an evaporatively-driven quasi-Roche lobeoverflow supersonic wind model for rotation-powered states. The mass loss must be low enough to escape IR/opticaldetection. The scenario is somewhat fine-tuned such that the companion wind is fast enough to inhibit a disk, whiledense enough such that angular momentum losses owing to turbulence are sufficient for gravitational influences tooverpower the pulsar wind. Such turbulence may also be driven by the radio absorption that predicates the eclipsingmechanism. Accordingly, R / R ∼ ( ˙ m g c / ˙ E SD ) (cid:29) R = 2 GM MSP /c is the Schwarzschild radius of the MSPand ˙ m g is the gravitationally-captured wind’s mass rate. However there are stability concerns – the pulsar terminationshock that arrests accretion flow and shrouds the pulsar and delineates the eclipsing medium may only be pushed outto a modest 2 −
100 multiples of pulsar light cylinder radius R LC (cid:28) a (Ek¸s˙I & Alpar 2005; Linares 2014) unless afeedback mechanism is operating. We show in an upcoming paper that if there is a feedback mechanism operating inRBs for the mass-loss rate from the companion, then such autoregulation may permit the shock to be stable muchfarther from the light cylinder out to orbital length scales for rotation-powered disk-free states (Wadiasingh 2017, inprep). A detailed discussion of the poorly-understood nuances of the irradiated companion mass loss, stability, andeclipse 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 afunction of inclination. This requires an a priori model for the shock geometry. For cases where the shock is orientatedaround the companion, two formalisms have been invoked for the free electron density underpinning the eclipses: anoptically 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 companionwind (Rasio et al. 1989, 1991; Tavani & Brookshaw 1991). We adopt the optically-thick formalism, as this leads toconstraints that are upper limits on R for a given i which we apply to B1957+20 as a test case in § θ max , R for agiven radio band. This is conceptually similar to radius-to-frequency mapping used in the study of pulsed emission inradio pulsars (Komesaroff 1970; Cordes 1978). Small values of θ max , R , corresponding to high radio frequencies, sampleregions 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). Instead we wish to constrain geometric shock parametersusing radio eclipses as a function of binary inclination i . This not only is of some utility to synthetic X-ray light curvesthat follow in §
4, but also informs on the ratio of wind ram pressures. Moreover, although the geometric model wepresent is somewhat degenerate on the parameter θ max , R , this motivates more systematic radio eclipse populationstudies of BWs and RBs.The general problem of an arbitrarily-shaped region occulting a source involves computational geometry techniquesthat may require inefficient ray casting or tessellating grids. To make the problem more analytically expedient, weassume an azimuthally symmetric form for the intrabinary shock with radial function { R ( θ ) | θ ∈ (0 , θ ∞ ) } along the axisof symmetry of the bow-shaped shock, with R ( θ = 0) = R , and generalize this approach to approximate the swept-backtails due to orbital motion. Azimuthal symmetry of the shock is expected to be an acceptable approximation for theintrabinary shock locale in the vicinity of the shock stagnation point, although overall the shock angle may be skewedrelative to the line joining the two stars yielding the extended-egress delayed-ingress radio eclipse phenomenology. Thisapproximation to the shock structure should especially be good when restricted to the ingress portion of eclipses aboutSC that are sharp and well-defined, unlike those that are more diffuse at egress and that are likely contaminated by0 Wadiasingh et al. plasma from the cometary-tail. In Appendix A, we develop an analytical formalism for eclipses of a point source byan arbitrary azimuthally-symmetric surface orbiting the source around a common barycenter, as well as parameterizethe 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 thatassume nonrelativistic, momentum-dominated winds neglecting gravity, through the balance of ram pressures, can befound in Wilkin (1996) and Canto et al. (1996). Here the companion shock, contact discontinuity, and shocked/deflectedpulsar wind are roughly spatially coincident compared with the length scale a , although not necessarily well-mixed, andinternal pressure contributions to the momentum-flux tensor are neglected. Internal pressure contributions generallyincrease the geometric thickness of the shock, particularly near the stagnation point where they are non-negligible,a complication that is largely peripheral to this Section. If the flow is not highly supersonic, the analytic forms aresignificantly 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-symmetricshocks can also be found in Gayley (2009) or Morlino et al. (2015), and generically result in increasing the shockopening angle.Since the physics and geometry of the pulsar wind and induced companion wind or magnetosphere are largelyunknown, we utilize a generic two-isotropic-colliding-winds analytic solution Canto et al. (1996) for the intrabinaryshock geometry that surveys varied shock asymptotic angles through a simple parameterization. This is readilyapparent from the gray contours in Figure 1 and qualitatively resembles hydrodynamic simulations of Tavani &Brookshaw (1991) for B1957+20. As explored in § R (cid:28) R iso ( θ ) = sin θ csc( θ + θ ) with θ cot θ = 1 + η w ( θ cot θ − , (6)and θ ∈ (0 , θ ∞ ) the polar angle defining the shock, with zero taken as the line separating the two stars and θ ∞ theasymptotic shock angle, θ ∞ − tan θ ∞ = π − η w , (7)where η w is the ratio of the two wind ram pressures and θ is an implicitly defined function of θ and η w . Explicitly, R is related to η w by R = √ η w √ η w . (8)We caution that the physical interpretation of η w may be misleading, especially for scenarios where the shock wrapsaround the pulsar where gravitational influences and unknown wind anisotropies are salient. Otherwise, when theshock is orientated around the companion, as for B1957+20 in § η w may be connected to the pulsar’s parametersif the wind of the companion is induced, e.g., Eqs. (10)–(11) in Harding & Gaisser (1990) and Eq. (9) in the nextSection.In §
4, we only employ the geometry of these analytical shocks rather than their physical velocity and densityprofiles, 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 fartherdown, while tangential velocity increases approximately linearly with the shock polar angle parameter θ . These arereadily apparent from a Taylor series expansion of analytical expressions in the thin-shell limit. At the stagnationpoint, the tangential velocity approaches a small value, if internal pressure contributions are small – this is indeedthe behavior observed in nonrelativistic hydrodynamic simulations (e.g., Figure 3 of Reitberger et al. 2014). Thetangential velocity v s and surface density Σ e variation with θ for hydrodynamic bow shocks can be shown to follow v s ∝ θ and Σ e = Σ (1 + wθ ) for θ (cid:46) w <
0, where | w | (cid:28) § hocks in MSP Binaries Figure 2.
Panels of fixed θ max , R depicting computed curves for the one-to-one coupling between eclipse fraction f E andthe shock stagnation point R for PSR B1957+20 in the axisymmetric optically-thick shock scenario with a range of orbitalinclinations i , as developed analytically in Appendix A. The panels depict successively higher values of θ max , R from left to right,with the two rightmost panels prescribing values of θ max , R that depend on R as stated. The bold blue curve highlights thebest-fit optical light curve modeling observational result i = 65 ◦ ± ◦ in Reynolds et al. (2007), while other curves survey theviable lower and upper systematic bounds as presented in van Kerkwijk et al. (2011). Note that ζ ≈ i ≈ ◦ is also found fromfavored outer-gap γ -ray pulsation fitting in Johnson et al. (2014). The mass ratio is fixed to q = 69 . R indicate the companion and volumetric Roche lobe R vL ( q ) radii (Eggleton 1983). Radio observationsfrom Ryba & Taylor (1991) are illustrated with the horizontal lines. The red region below 600 MHz, where eclipses are to somedegree symmetric, isolates roughly the region of validity of the axisymmetric calculation. density profile is probably an inaccurate proxy for the spatial distribution of particle acceleration and cooling; thus amore general procedure is also described in § 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 geometryin Eq. (6) using Eq. (A15). This is conditional on the crucial parameter θ max , R , where the medium transitions fromoptically thick-to-thin for a given observing frequency. Note that R and i are the geometric parameters that areindependent of observer frequency, therefore θ max , R is the only parameter in the model that connects to the frequencydependence of symmetric eclipses.In Figure 2, we display the eclipse fraction f E dependence on R of PSR B1957+20 with a fixed mass ratio q = 69 . ζ ≈ i ≈ ◦ found from Johnson et al. (2014). The axisymmetric computations for constraining R should be more accurate for the eclipses at ν (cid:38)
600 MHz that are largely symmetric about SC; this is indicatedby the red region in the panels. We neglect the minor degenerate observational coupling between mass ratio andinclination angle due to uncertainties in the optical mass function. The four panels depict successively larger values ofthe parameter θ max , R left to right. In the leftmost two panels the value of θ max , R = { π/ , π/ } is independent of R ,while the rightmost two panels impose values that depend on R , e.g. through Eq. (7). In particular, the prescription R ( θ ) sin θ = 4 R selects the value of θ max , R for a given R such that the transverse shock length scale is 4 R . Onsuch a scale far from the shock head, hydrodynamic instabilities will likely develop that may break the assumptionsin our rudimentary model. The rightmost panel chooses an exceptionally large value of θ max , R = 0 . θ ∞ that servesas an extreme limit for what θ max , R may be and represents a very substantial occluding volume. Such large values of θ max , R > π/ θ max , R , these values are included for completeness. Similarly, values much smaller than θ max , R < π/ f E tending to R → . L point, also a rather unreasonable scenario that requires η w ∼ § f E , it is clear from Figure 2 that larger values of θ max , R are compatible with smaller values of R . For large inclinations near edge-on, there are clearly constraints on θ max , R ; in particular for i = 85 ◦ , θ max , R (cid:46) π/ R > R ∗ . On the other hand, for i = 65 ◦ , the constraint is looser with θ max , R > { π/ , π/ , R ( θ ) sin θ = 4 R , . θ ∞ } Wadiasingh et al.
Figure 3.
Schematic representations of the PSR B1957+20 system, to scale, projected on the plane of the sky for inclinations55 ◦ − ◦ (columns), q = 69 .
2, and orbital phase 0.034 from SC at eclipse ingress/egress. Rows of insets indicate several fixed θ max , R , the same values as those in Figure 2, with R chosen such that f E ≈ R ∗ = 0 . R vL ≈ . R < R ∗ is obligatory are omitted. For all panels,the color coding emphasizes, schematically, the locales of first-order Doppler boosting. A velocity profile of v s ∝ θ tangent tothe shock is imposed; the blue or red coloring accentuate those regions of the shock where the v s along the shock is toward oraway the observer line-of-sight, respectively, with intensity of coloring scaled with projected velocity component magnitude fromEq. (25). This coloring qualitatively demonstrates the geometrical contribution to the emissivity integral in § corresponding to limits R (cid:46) { . , . , . , . } , respectively, for f E (cid:46) R are modestlymore stringent for ν >
600 MHz for the same range of θ max , R . Therefore we conclude that for i = 65 ◦ , R ≈ . − . R ≈ .
235 for θ max , R = π/ R ≈ . § R leadsto DP light curves that have too-wide of a peak separation as will become apparent in due course. Some geometricrealizations for f E ≈ ν ≈
600 MHz, are illustrated in Figure 3 with corresponding numericalvalues of θ max , R = { π/ , π/ , R ( θ ) sin θ = 4 R , . θ ∞ } (columns) and R . Some geometric trends are clearly evidentin Figure 3, for instance larger i requiring smaller R for the same f E . Similarly, larger θ R , max for fixed i allows forlower R .The upper limits found for R allow for an estimate of the companion wind pressure if ˙ E SD is known. We mayexpress η w , with P the wind pressure due the companion, as η w = P δ Ω ˙ E SD / (4 πc ) ∼ . π ) R ∗ σ B T /c + ξ ˙ E SD /cδ Ω ˙ E SD / (4 πc ) (cid:28) , (9) hocks in MSP Binaries T cold is the isotropic unheated temperature of the stellar companion (i.e., the intrinsic unirradiated radiationpressure from the secondary), δ Ω / (4 π ) is the fractional isotropic pulsar wind solid angle subtended by the shock. Ifa strong magnetosphere is not the principal source of ram pressure against the pulsar wind, then the parameter ξ embodies the fractional energetic efficiency, in units of ˙ E SD , of the induced companion wind generated by unspecifiedprocesses. The R found above are compatible with the thermally-driven wind or companion magnetosphere scenariosof Harding & Gaisser (1990). Numerically, for typical values R = 0 . − . η w is between 6 to 18% by inverting Eq. (8). From Eq. (9), we can form an estimate of the energeticefficiency ξ of the induced wind, if the intrinsic or induced magnetic field of the companion is small, ξ = 14 π (cid:18) R − R (cid:19) δ Ω − πR ∗ σ B T ˙ E SD ≈ π (cid:18) R − R (cid:19) δ Ω , (10)where the ratio of cold intrinsic stellar to pulsar power can be neglected in BWs and RBs, since it of the order ∼ − − − , and thus ξ is a simple function of stagnation point R . This equation is unphysical for R → . θ ∈ (0 , π/
2) which mayparticipate in the heating of the companion can be routinely found δ Ω / (4 π ) ≈ R (1 + 2 R ) / R (cid:28) ξ is of the order 0 . − R . This efficiency is similar in order-of-magnitude to the hemispherical quiescent induced photosphericheating fraction, i.e. 2 πR ∗ σ B T / ˙ E SD ∼ .
4% for B1957+20 where T hot ≈ B ∗ may be found by assuming companionpressure is entirely due to a magnetosphere at the stagnation point (e.g., Eq. (23) of Harding & Gaisser 1990), B ∗ (cid:46) × (cid:32) ˙ E SD erg s − (cid:33) / (cid:18) cm a (cid:19) (cid:18) . R ∗ /a (cid:19) (cid:18) R . (cid:19) (cid:18) . − R (cid:19) G , (11)which is relatively small compared to typical surface fields found in degenerate cores, and comparable to kilogauss fieldsfound in T Tauri stars (Johns-Krull 2007). In principle, such fields may be detectable in future, but currently requiresrelatively bright (magnitude (cid:46) −
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 fieldconstraint would be useful in constraining the physics of the pulsar and induced companion winds, as well as to appraisethe Applegate (1992) model for a tidally-driven convective dynamo. Indeed, a strong companion magnetosphere wouldalso dramatically alter shock MHD jump conditions, impact the relevant acceleration mechanisms in the shock, andinfluence 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 ◦ and85 ◦ , with the former only sampling regions of the confined flow peripheral to the shock head. This segues to theinclination-dependent numerical computation, using Eq (A15) of the growth rate of the eclipse fraction f E as functionof θ max , R , depicted in Figure 4. The appropriate interpretation of Figure 4 should be restricted to the symmetricpart of eclipses, i.e. below about f E (cid:46) .
1. The curves terminate at θ ∞ , and we note the computed f E < . θ max , R → θ ∞ . For high inclinations,where the eclipses largely sample the shock head, the growth rate is approximately linear with θ max , R contrasting thenonlinear growth rates for lower inclinations which sample the periphery and tail of the shock. If the optical depth τ due to scattering or absorption by a cross section σ ν ∝ ν − m is given by 1 (cid:28) τ ∝ σ ν Σ e , R for column density Σ e , R ,then the spatial variation of the column density can be probed. Given growth curves f E ≡ g ( θ ), for an observedfrequency dependence in radio eclipses f E ∼ ν ( θ ) − n , as observed for B1957+20 and RB J2215+5135 with f E ∝ ν − . ,then ν ( θ ) ∼ g − /n and the spatial distribution of the integrated column for a given optical depth is proportional toΣ e , R ( θ ) ∼ g ( θ ) − m/n . For instance, if g ( θ ) ∼ θ l , n = 0 . m = 2 for free-free absorption in the Rayleigh-Jeans limit,then d log Σ e , R ( θ ) /d log θ ∼ − l , an inclination-dependent line-of-sight column density that can attain a rather steepand sharp profile. This motivates frequency-dependent radio population studies of similarly eclipsing BWs and RBsfor inclination-dependent trends of f E ( ν ). This is of generally low utility to the X-ray Doppler-boosted emission in § f E growth curves for B1957+20 and other SC-centered spiders. Such a parallel-wind geometry is self-4 Wadiasingh et al.
Figure 4.
Growth of f E as a function of θ max , R ∈ { , θ ∞ } illustrated for six values of R between 0 . − .
35 at fixed i forthe scenario where the shock is around the companion. The curves terminate at θ ∞ , resulting in a maximum f E < . R and i . The flat portions at low θ max , R for i = 85 ◦ and 90 ◦ are due to the companion of R ∗ ≈ . similar in the R (cid:28) R has no impact on theshock opening angle but only serves as a scaling parameter. Since the companion and pulsar winds are not anticipatedto be isotropic but somewhat stronger about the line joining the two stars in the orbital plane, this geometry is arealization where wind anisotropies obligate a much smaller shock opening angle. Such wind anisotropy is expected foran anisotropically-irradiated companion whose evaporatively-driven wind is only influential on the day-side hemisphereof the companion. For this auxiliary geometry, it is found that the constraints on R are systematically larger thanthe isotropic-winds geometry, and less sensitive to the value of θ max , R since f E plateaus in the far-downstream tailof the shock. Since the shock angle is small, an additional upgrade to ballistic tail sweepback and eclipse asymmetryis also explored. This parameterization of eclipse asymmetry in terms of outward wind speed serves as an additionalmotivation 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 ( θ ) sin θ is amonotonically rising function of θ , as is the case for all shocks in the limit of small sweepback, the eclipse fraction isindependent of shock geometry for any given binary inclination and only depends on the shock’s largest attained polarangle θ max , R . For cases where the shock is swept back due to Coriolis effects which may introduce φ dependences andbreak the axisymmetry of the shock, the axisymmetric estimate below is a lower limit on the shrouding fraction (seeAppendix A). For the eclipsing medium, the only assumption here is that it lies beyond the MSP termination shock,i.e. the geometric complement of the cavity excavated by the pulsar wind. Therefore, shrouding at IC is suppressedand only anticipated at inclinations away from edge-on in this model. The stand-off distance R is measured from thepulsar, rather than the companion, in this context.If the shock is at an angle ϑ sb with respect to the orbital angular momentum vector, the eclipse fraction f E is relatedto θ max by routine spherical trigonometry,cos θ max , R = cos i cos ϑ sb + sin i sin ϑ sb cos( πf E ) . (12) hocks in MSP Binaries i = ° ° ° ° ° ° θ ∞ | R , m a x = . θ ∞ | R , m a x = . θ ∞ | R , m a x = . θ ∞ | R , m a x = . π π π π π π θ max, R f E Figure 5.
Analytical eclipse fraction f E growth curves for any DP IC-centered system where the shock surrounds the pulsar,given by Eq. (13), are shown highlighting the constrained inclination of J1023+0038 in the green to purple curves. The redhorizontal lines illustrate typically observed asymmetric total f E from Archibald et al. (2009, 2013). The dark yellow verticallines show the maximum R that may accommodate the asymptotic shock opening angle θ ∞ for the specific geometry of Eq (6). We adopt the case ϑ sb = π/ ϑ sb = π/ ζ ≈ i and the absence of jet-like anisotropy in the MSP wind (e.g., Coroniti 1990;Lyubarskii 1996; Bogovalov 1999), a reasonable conjecture at this juncture. This yields f E = 1 π arccos (cid:18) cos θ max , R sin i (cid:19) (13)for π/ − i ≤ θ max , R ≤ π/ i with the lower and upper limits corresponding to f E = 0 and 1 respectively, in contrast toFigure 4 which can never exceed 0 .
5. Thus radio eclipses in MSP binaries where X-ray emission is IC-centered constrainthe shrouding by the maximum shock polar angle θ max , R in a model-independent manner, and values of the polar angleare associated with a physical integrated line-of-sight density profile. Note that this formalism allows for axisymmetricshocks that are skewed at an angle ∆ φ s relative to the line joining the two stars. If φ in and φ eg are the orbital phasesof ingress and egress eclipses relative to SC, respectively, then the shock asymmetry is directly measurable at a givenfrequency (corresponding to a θ max , R in this model) as the antisymmetric part (1 / φ eg − φ in ) ≈ ∆ φ s / (2 π ).The contours corresponding to Eq. (13) are illustrated in Figure 5 with approximate J1023+0038 observational radioeclipse fractions from Archibald et al. (2009, 2013) for a range of inclinations. For parallel-wind shocks, large valuesof θ max , R > π/ L z (cid:29) R where hydrodynamic instabilities are expected tobe influential. Moreover, as we show in § R corresponding to θ ∞ are shown in Figure 5, thus, one may constrain the upperlimit of the shock opening angle based on the largest stable eclipse fractions observed. For J1023+0038, from eclipsesat 150 MHz of f E, (cid:38) . R (cid:46) .
4. We caution that in such IC-centered spiders, the naivephysical interpretation of R constraints in terms of wind ram pressures given by Eqs. (6)–(8) is erroneous due thecommanding gravitational influence of the MSP past the L point, and depends on the specifics of angular momentum6 Wadiasingh et al. loss of the companion baryonic wind. Instead, R may be envisaged as a convenient parameterization of the shockopening angle. THE DOWNSTREAM BULK LORENTZ FACTOR Γ, SHOCK MIXING, AND BARYON LOADINGFor non-relativistic shocks in Eq. (6) and Eq. (B17), the tangential velocity flow increases approximately linearlywith the shock polar angle θ or symmetry axis z for points close to the stagnation point where the shock is very nearlya spherical cap. Lacking a self-consistent relativistic MHD shock geometry, for simplicity we extend this flow scalingrelation to the mildly relativistic regime assuming a scaling for the bulk specific relativistic momentum, p Γ ≡ Γ β Γ = (Γ β ) max (cid:18) θθ max , X (cid:19) , (14)with θ max , X corresponding to a characteristic angular scale where the shocked pulsar wind is pressure-confined inthis manner and defines the region of interest where particle acceleration and synchrotron processes operate. In thisrudimentary model the choice of θ max , X , which is independent of the radio counterpart θ max , R , is set to a benchmarkvalue π/ § β Γ is simply the spatial part of the4-velocity that appears in the pressureless relativistic fluid continuity equation or energy density tensor. The fluidspeed and Lorentz factors along the shock are a function of θ , viz. β ( θ ) = (Γ β ) max θ (cid:113) (Γ β ) θ + θ , X , Γ( θ ) = 1 (cid:112) − β ( θ ) . (15)We adapt this ad hoc spatially-dependent flow speed prescription Eq. (14)–(15) for the calculationand parameterizationof Doppler factors to synthesize light curves in § c/ √
3. The physicalinterpretation of the bulk Lorentz factors Γ depends on the baryon loading of the shocked pulsar wind by the shockedcompanion component across the contact discontinuity, and couples to the companion mass loss rate in a nontrivialmanner; a loading n e /n i ∼ is sufficient for ions to be influential where n e,i are electron and ion number densities. Forthe energy budget ˙ E SD , “well-mixed” ion-dominated hypothesis such as that of Bednarek (2014) requires a substantiallylower companion mass loss rate than their assumed ˙ m ≈ . × g s − to yield bulk Lorentz factors large enoughto yield DP X-ray modulation as observed, if attributed to Doppler boosting in B1957+20 (cf. § E SD since it augmentsthe leptonic population in the shocked pulsar wind, and is consequential for spectroscopically resolved light curvesin § § β ) max and otherparameters in a fixed energy bin. Notice that the level of mixing will also influence the geometric thickness acrossthe contact discontinuity of the leptonic pulsar and baryonic companion winds. This is because in a kinetic pictureshock thickness is proportional to the mean Larmor radius of the particles, while in a macroscopic fluid picture, thegeometric thickness of the shock structure is related to the internal pressure terms which are anticipated to be largerfor the nonrelativistic baryonic component versus the momentum-dominated pulsar leptonic component. Thereforeone may expect mixing to increase the geometric thickness and volume where leptons radiate.We anticipate a partitioned shock morphology in spiders, with a shocked companion wind and pulsar wind, separatedby a contact discontinuity. From flux-freezing arguments in an idealized steady-state laminar MHD framework, thereis no charged particle transport across the contact discontinuity except perhaps at the stagnation point. This lackof transport seems to be in tension with the assumptions of Romani & Sanchez (2016). RMHD instabilities andkinetic-scale processes alter this overly-simplified picture for mixing. Following Morlino et al. (2015), an estimate ofmixing due to kinetic-scale effects may be discerned by consideration of cross-field diffusion of charges. In the Bohmlimit of strong turbulence and wave-particle interactions, the perpendicular diffusion coefficient κ ⊥ proceeds at thegyroscale attaining κ ⊥ (cid:46) κ (cid:107) = r L v/
3. Although such strong turbulence is not anticipated in the shocked companionwind, this estimate constitutes a lower limit for the cross-field diffusion timescale τ D , ⊥ of protons, τ D , ⊥ (cid:38) (cid:96) κ (cid:107) = 1 . × (cid:18) (cid:96) cm (cid:19) (cid:18) B s (cid:19) (cid:16) v p cm s − (cid:17) − s , (16) hocks in MSP Binaries v p = 10 cm s − (cid:29) a/P b for the proton speed and assumed (cid:96) = 0 . R ∼ cm, the lengthscale corresponding to the separation between the shock components, roughly that recovered from RMHD simulations(Bucciantini et al. 2005). The value of τ D , ⊥ attained in Eq. (16) is much greater that the typical advection timescaleof ions of 10 −
100 s in the shocked companion wind, thus mass loading due to diffusion is unlikely to be a significantinfluence unless σ is significantly smaller and the wind is particle dominated. Eq. (16) adapted to pairs in the shockedpulsar wind is a factor of 1836 larger, thus less constraining. The onset of significant baryon loading due to kineticeffects then corresponds to (cid:104) B s (cid:105) (cid:46) − and σ (cid:46) − in the Kennel & Coroniti (1984) prescription for obliqueshocks, values that are largely excluded since they obligate a Lorentz factor of leptons in the shocked pulsar wind γ e (cid:38) (which violate the Hillas bound for containment in the shock) to emit at the critical synchrotron frequencycorresponding to the >
10 keV NuSTAR band. The contribution of mixing due to instabilities is difficult to quantifyand depend on σ in a nontrivial manner, but less germane for σ (cid:28) R .Additionally, using the magnetic field constraint Eq. (5) and the pair wind density Eq. (3), it is simple to form anupper limit for the shocked magnetization parameter σ shocked if there is no baryon loading, σ shocked (cid:46) B , min πn e , MSP (cid:104) γ e (cid:105) m e c ≈ × M ± (cid:104) γ e (cid:105) (cid:18) B s , min (cid:19) (cid:16) r s cm (cid:17) (cid:18) erg s − ˙ E SD (cid:19) / . (17)Therefore, modest values of M ± and (cid:104) γ e (cid:105) are sufficient to yield σ shocked (cid:46) R . Therefore, the mixing of shock components is insignificantunless small-scale instabilities are prolific near the head of the shock or there are unaccounted for magnetic fieldsquasi-parallel to the shock normal. Mixing with the companion wind electron/ion plasma also dilutes the energybudget for the shock-accelerated energetic positrons that escape the system. Hence, models that posit MSP binaryshocks as significant astrophysical sources of energetic cosmic-ray positrons at Earth assume low mixing (e.g., Venteret al. 2015). Unless the unknown mass loss rate of the companion is small, low levels of mixing are also necessaryfor the mildly-relativistic bulk Lorentz factors ascribed by Eq. (14)–(15) and in § § GEOMETRIC DOPPLER-BOOSTED ORBITALLY MODULATED SYNCHROTRON EMISSION4.1.
Formalism
We assume local quasi-isotropy in pitch angle of the differential electron number density distribution along the shockin the comoving frame of the bulk flow in a thin emitting region such that ¯ n e ( ¯ r , γ e ) → Σ e ( θ, φ, γ e ) / ∆( θ, φ ) at eachlab-frame coordinate pair { θ, φ } on the shock, where ∆( θ, φ ) is the shock thickness and Σ e the differential electronsurface density distribution. In reality, the local electron distribution accelerated in relativistic shocks by DSA orprolific magnetic reconnection is expected to be highly anisotropic relative to the mean field direction in the comovingframe. However, for locales near the stagnation point where the residence timescale is large compared to the adiabaticand convection timescales, the flow is expected to be well-isotropized with the field randomly oriented, especiallyif magnetic dissipation is the principal particle acceleration mechanism. Moreover, for locales where the convectionspeed 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 theunderlying electron pitch angle distribution.The differential volumetric emissivity, in the comoving frame (denoted by bars over variables), is defined by¯ j (cid:15) f = d ¯ Ed ¯ V d ¯ td ¯ (cid:15) f d ¯Ω f = m e c ¯ (cid:15) f ¯˙ n isoSR , (18)where ¯˙ n isoSR is the isotropic differential synchrotron photon production rate at each interaction point. The variables V , t , (cid:15) f and Ω f constitute the volume, time, outgoing photon energy, and outgoing photon solid angle, respectively, withphoton energies in units of m e c .8 Wadiasingh et al.
For a power-law distribution of electrons with index p between Lorentz factors γ min and γ max and spatial range θ ∈ (0 , θ max , X , denoted by the Heaviside step function Θ, such that ¯ n e ( γ e , θ ) = Σ e ( θ, φ ) / ∆( θ, φ )( γ e /γ min ) − p × Θ( γ e γ min , γ max )Θ( θ ; 0 , θ max , X ), the emissivity in the synchrotron power-law regime γ min (cid:28) (cid:112) ¯ (cid:15) f B cr / ¯ B s (cid:28) γ max , isgiven by Dermer & Menon (2009),¯ j ¯ (cid:15) f ≈ C ( p ) (cid:18) (cid:19) ( p +1) / Σ e ( θ, φ )∆( θ, φ ) σ T c U cr (cid:18) ¯ B s ( θ, φ ) B cr (cid:19) ( p +1) / ¯ (cid:15) − ( p − / f , (19)where ¯ B s ( θ, φ ) is the post-shock turbulent magnetic field magnitude in the comoving frame, U cr = B / (8 π ), B cr ≈ . × G the quantum critical or Schwinger field, and C ( p ) = 2 ( p − / √ √ π (cid:20) Γ (cid:18) p (cid:19)(cid:21) − Γ (cid:18) p − (cid:19) Γ (cid:18) p (cid:19) Γ (cid:18) p (cid:19) . (20)In the idealized MHD formulation where the bulk flow is force-free and there is no magnetic field perpendicular tothe flow, we note that Lorentz transforming from the lab frame to the comoving blob frame preserves the parallelcomponent (cid:104) B s (cid:105) .To suitably develop light curves that would correspond to an observable, we form the differential synchrotronluminosity, L SR (Ω f , (cid:15) f ) ≡ d L SR / ( d Ω f d(cid:15) f ). This is the total source luminosity L SR , or volume-integrated emissivity,binned in solid angle and energy elements. We take advantage of the simple Lorentz transformation property of thespectral emissivity ¯ j (cid:15) f , namely that j (cid:15) f /(cid:15) f is invariant under boosts. Forming the differential luminosity in the observerframe, 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: L SR (Ω f , (cid:15) f ) ≡ d L SR d(cid:15) f d Ω f = (cid:90) dV (cid:15) f (cid:32) j (cid:15) f (cid:15) f (cid:33) = (cid:90) dV δ ¯ j (cid:15) f , (21)with the emissivity and photon energies being computed using the Doppler-shifted photon energy¯ (cid:15) f = (cid:15) f /δ D . (22)Observe that the factor δ in Eq. (21), rather than δ found routinely in relativistic jet contexts, arises becausethe integration over the emitting volume element is chosen to be in the observer’s frame rather than the comovingframe. This orbitally-modulated Doppler factor δ D is calculated for each point along the shock and dependent on theprescribed local bulk speed β Γ ( θ ) and bulk Lorentz factor Γ( θ ) along the shock, defined by δ D ( θ, φ, Ω b t, i ) = 1Γ( θ ) (cid:0) − β Γ ( θ ) ˆ n v · ˆ u (cid:48) (cid:1) , (23)where ˆ u (cid:48) is the unit tangent vector along the polar direction of the shock in the inclined and orbital phase-rotatedcoordinate system (cf. Appendix A for definitions and conventions). Thus, we have u ≡ dR x dθ ˆ x + dR y dθ ˆ y + dR z dθ ˆ z , ˆ u = u u (24) ˆ u (cid:48) = Λ i Λ Ω b t ˆ u . (25)In addition, ˆ n v is the unit vector in the direction of the observer, applicable to all orbital phases. Throughout, weemploy the scaling Eq. (14)–(15), and associated β 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 iswhat defines the relative sharpness of peaks in the ensuing figures of orbital modulations of the flux. For the assumedsymmetry, the integration Eq. (21) is two-dimensional in variables θ and φ 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 dV ≈ ∆( θ, φ ) dA = a ∆( θ, φ ) R ( θ ) sin θ (cid:112) R ( θ ) + ( dR/dθ ) dθ dφ. (26) hocks in MSP Binaries
19. The factor ∆( θ, φ ) cancels with the factor in Eq. (19), and the net result is that the orbital-modulated differentialluminosity or intensity is proportional to L SR ∝ (cid:15) − ( p − / f (cid:90) dθ dφ R ( θ ) sin θ (cid:112) R ( θ ) + ( dR/dθ ) δ p − / ς ( θ, φ ) (27) ς ( θ, φ ) = Σ e ( θ, φ ) (cid:2) ¯ B s ( θ, φ ) (cid:3) ( p +1) / (28)when p is a constant along the region of interest. Crucially, the flux ratio of Eq. (27) generates energy-independent lightcurves when p is spatially independent. In this regime, the ensuing large curves are largely regulated by geometricinfluences of the observer and ascribed velocity profile Eq. (14) through Eq. (23). It is quickly discerned that themaximum amplitude of orbital modulation is attained when ˆ n v · ˆ u (cid:48) ≈ /β Γ realizing a flux enhancement of order (cid:46) p − / modulo weighting of the emissivity across the full shell and observer impact angles with respect to thetangential velocity component in Eq. (27) that moderates this upper limit.The accelerated/cooled electron power-law and magnetic field spatial distribution in the thin shell, encapsulated in ς , is poorly understood. Geometric and Doppler boosting influences must dominate in the integral Eq. (27) far fromthe shock nose since this is a necessary condition for the existence and generation of DP light curves. Therefore theweighting ς is expected to be a modest nuance on the light curve morphology, pending a future proper treatment ofself-consistent particle transport. Even so, the underlying particle acceleration mechanism is unknown which couples tothe form of ς in a transport model, and the unknown spatial distribution of p . Crucially, if p is spatially dependent, theoverall spectral index may deviate from the usual ( p − / § ∂ς/∂θ = ∂ς/∂φ = 0 with p = 2 and routinely integrate Eq. (27)for the isotropic-winds geometry of Eq. (6) for a single thin shell in the IC-centered scenario, IC corresponding toangular phase π . For such a shock orientation, geometric occlusion by the companion is a negligible influence. A phaseplots atlas of light curves, depicted in Figure 6 normalized to unity at SC, establishes several features as a functionof parameters R , (Γ β ) max , and i , with θ max , X = π/
2. Observe that with a single shell, only a single broad peak orDPs may be produced, as discernible in the “U” or “V” patterns in the top two rows. The boost parameter (Γ β ) max straightforwardly regulates the relative sharpness and width of the peaks. The characteristic shock asymptotic openingangle, parameterized by R , regulates peak separation at fixed i with larger shock opening angles (corresponding tolarge R ) yielding wider peaks. The characteristic shock opening angle is also regulated by the weighting ς . Transientchanges in X-ray peak separation therefore may be interpreted as changes in the effective shock opening angle withrespect to the observer or an equivalent spatial change in ς , since emission near the shock nose is not Doppler-boosted.There is degeneracy in peak separation with i and R highlighting the need for multiwavelength constraints on i . Theother degeneracy for the peak width between R and (Γ β ) max at fixed i may be resolved in a future particle transportand mixing model, since (Γ β ) max cannot be arbitrarily large by energy budget considerations. If there is universalityamong BWs and RBs in the shock structure or opening angle, correlations between X-ray peak separation and i maybecome evident across a population of BWs or RBs.Note that the assumption of pitch angle isotropy also suppresses a phase-dependent polarization signature. Apolarization-dependent calculation follows routinely from the development for the synchrotron emissivity, but requiresa specification of the unknown magnetic structure and level of turbulence. Future proposed X-ray polarimetry instru-ments (PolSTAR: Krawczynski et al. 2016; XIPE: Soffitta et al. 2013; IXPE: Weisskopf et al. 2016) may be able todiscriminate among types of acceleration mechanisms from orbital phase-dependent Stokes parameters. High levels oflinear polarization, for example, may suggest radiative losses in an ordered magnetic field rather than turbulent recon-nection, with orbital phase-dependent polarization yielding tomographic information of the magnetic field geometryon the shock. 4.2. Application to PSR B1957+20, an SC-centered System
Using Eq. (6) we compute the volume-integrated emissivity ratio of superior-to-inferior conjunction in Figure 7,and take ς constant or ς ( θ ) = ς (cid:2) − ( θ/θ max , X ) (cid:3) at a fixed energy (cid:15) f where the particle power-law is valid, with allnumerical constant factors canceling. We choose p = 2 as a benchmark value, which is somewhat harder than impliedfrom photon indices Γ X ∼ . − ς ( θ ) is known from hydrodynamic bow shocks, but likely does not trace the true time-averaged emitting particledistribution. The integration is taken to θ max , X = π/
2, that is, we only consider the head of the bow shock that0
Wadiasingh et al.
Figure 6.
Partial atlas of IC-centered light curves for synchrotron emission, normalized to unity at SC (white coloring),computed for p = 2, θ max , X = π/ ς = constant using Eq. (27). R is measured from the MSP. The top, middle and bottomrows depict the i , R and (Γ β ) max dependency on the synthetic light curves, respectively. The columns contrast four pairingsof parameters somewhat extreme in range. No shadowing by a companion is included. participates for the hemisphere of the companion facing the pulsar. This portion of the shock geometry is expected tobe largely axisymmetric and comprise the preponderance of accelerated charges and emission.We employ the results of § R for a given inclination in the axisymmetric case. Forlower maximum bulk speeds β max = 0 .
5, the left panel in Figure 7, although there is orbital modulation it is relativelymodest and flat around SC. Little to no DP structure is exhibited except when the shock stagnation point is taken atthe companion surface, assumed to be 90% of the Roche Lobe radius for B1957+20, where the influence of geometricoccultation of the shock by the companion (solid curves) is larger. However, this DP morphology is rather flat andwould require a large fraction of emission to be concentrated near the nose to produce the observed peak-dip ratioin Huang et al. (2012), and even then would not produce the correct ≈ . R ≈ R ∗ would also appear to be in tension with the radio eclipse estimates in § R , the influence of shadowing (solid versus dashed curves) is smallespecially for the more moderate inclination i = 65 ◦ . In general, shadowing of unboosted optically-thin emission canonly produce a single dip, not DPs; high maximum bulk velocities, greater than c/ √
3, are required to establish theDP structure. Therefore, we conclude that shadowing of the emission region by the companion alone cannot explainthe DP light curve features, particularly peak separation and peak width. hocks in MSP Binaries % β max = i = o , R = ϛ ( θ ) i = o , R = ϛ ( θ ) i = o , R = ∂ϛ ( θ ) ∂θ = i = o , R = ∂ϛ ( θ ) ∂θ = i = o , R = R * , ϛ ( θ ) i = o , R = R * , ∂ϛ ( θ ) ∂θ = - π - π π π Angular Orbital Phase R e l a t i v e L u m i n o s i t y % β max = i = o , R = ϛ ( θ ) i = o , R = ϛ ( θ ) i = o , R = ∂ϛ ( θ ) ∂θ = i = o , R = ∂ϛ ( θ ) ∂θ = i = o , R = R * , ϛ ( θ ) i = o , R = R * , ∂ϛ ( θ ) ∂θ = - π - π π π Angular Orbital Phase R e l a t i v e L u m i n o s i t y Figure 7.
Orbitally-modulated synchrotron flux ratios of superior-to-inferior conjunction for β max = { . , . } with θ max , X = π/
2, at an arbitrary energy where the power-law approximation is valid, for different inclinations and shock stand-off R generallyconsistent with radio eclipses. Phase zero defines SC of the pulsar. The solid and dashed curves are cases, where the shock isshadowed and unshadowed by the companion , respectively. θ max,X = π / i = ° R = β max = β max = β max = N o r m a li ze d I n t e n s i t y i = ° R = i = ° R = i = ° R = i = ° R = i = ° R = N o r m a li ze d I n t e n s i t y i = ° R = i = ° R = i = ° R = i = ° R = i = ° R = π π π π Angular Orbital Phase N o r m a li ze d I n t e n s i t y i = ° R = π π π π Angular Orbital Phase i = ° R = π π π π Angular Orbital Phase i = ° R = π π π π Angular Orbital Phase i = ° R = π π π π Angular Orbital Phase
Figure 8.
Synthetic light curves highlighting the R and inclination dependence of the DP structure, with θ max , X = π/
2. Thefluxes are scaled to unity at maximum for clarity. The solid and dotted curves correspond to ς = constant and ς ( θ ) prescriptions,respectively. All light curves are computed for p = 1 . X = 1 . Larger bulk speeds in the right panel in Figure 7 do produce characteristically DP synthetic light curves, whichare in qualitative agreement with Huang et al. (2012) for the DP positions. As expected, modulation increases withinclination i for either shock geometry. Given the large error bars on the data in Huang et al. (2012) and the earlydevelopment stage of our model, we do not attempt any fitting. There is a complex interplay between the shockgeometry, the system inclination, the value of R (which principally moderates the influence of shadowing), the bulkLorentz factor at points along the shock, and the surface density profile ς ( θ ). Smaller inclinations reduce the peakseparation for a given set of parameters, which can also be mimicked by altering the prescribed shock geometry forwider opening angles. Shadowing is generally a negligible influence, but may modestly enhance the dip at SC for i = 85 ◦ .It is clear from this quantitative parameter exploration for B1957+20 that higher bulk Lorentz factors are preferredsuch that β ≈ . − . ς ( θ ) prescription, is essential tomitigate the larger than observed amplitude and peak-dip ratio.4.3. Application to PSR J1023+0038, an IC-centered System Wadiasingh et al. θ max,X = π / i = ° R = β max = β max = β max = N o r m a li ze d I n t e n s i t y i = ° R = i = ° R = N o r m a li ze d I n t e n s i t y i = ° R = i = ° R = π π π π Angular Orbital Phase N o r m a li ze d I n t e n s i t y i = ° R = π π π π Angular Orbital Phase
Figure 9.
Synthetic light curves highlighting the R and inclination dependence of the DP structure, with θ max , X = π/ ς = constant and ς ( θ ) prescriptions, respectively. All light curves are computed for p = 1 . X = 1 . The observed X-ray light curve of PSR J1023+0038 is centered at pulsar IC, so that the shock must be surroundingthe MSP instead of the companion. In this case, we take R = 0 is at the MSP rather than at the companion. In Figure 8we compute the symmetric Doppler-boosted light curves with p = 1 . θ max , X = π/
2, normalized to peak maximum,for a variety of inclinations and β max = { . , . , . } , the largest value corresponding to a maximum bulk Lorentzfactor of about 5. This calculation constitutes a variant of Figure 6, focusing on parameters relevant to J1023+0038restricted to i = 35 − ◦ , and p = 1 . X ≈ .
1. Moreover, the ς ( θ ) = ς (cid:2) − ( θ/θ max , X ) (cid:3) prescription is illustrated with dotted curves. Since the shock surrounds the MSP, no shadowing by the companion isincorporated in the calculation since even a fully Roche lobe-filling companion for the given low inclinations 35 ◦ − ◦ only obscures a negligible fraction of the shock surface emissivity, and then so only near SC for emission that isDoppler de-boosted. As for B1957+20, we do not attempt any fits here, pending the development of a more completeand self-consistent model of particle transport, cooling and shock geometry, but rather explore shock geometry andphysical parameters in a more generic way that will be useful to such a future study.A wide range of shock opening angles are surveyed by the range of R values. Smaller values of R , corresponding tonarrow bow shocks, for instance of the parallel-wind type, are generally disfavored by failing to yield DP modulationeven the largest prescribed value of β max = 0 .
98 except at the highest inclinations (corresponding to unusually lowpulsar masses) for either ς constant or ς ( θ ) = ς (cid:2) − ( θ/θ max , X ) (cid:3) . Larger values of β max than considered here mayallow these smaller values of R , a degeneracy in the model. The bulk Lorentz factor parameter largely controls thepeak width and the depth of the dip at IC. It is evident that R (cid:38) . § R (cid:46) .
4. For the intermediate inclination i = 45 ◦ corresponding to an MSP mass of M MSP ≈ . M (cid:12) employing the radial velocity K value from McConnellet al. (2015), we may constrain R (cid:46) . β max > .
8. Higher inclinations can yield DP structure with lower β max values. The variation in peak separation with R implies that decreasing or increasing peak separations may beevident in X-ray light curves of RBs promptly preceding or following an LMXB state transition, respectively.Interestingly, for large values of β max , the ς ( θ ) concentration near the stagnation point does enhance DP modulationfor small values of R , but has the opposite effect for larger values. This behavior for small R results from differentialweighting and contributions of the observer line-of-sight cutting across the beamed portion of the shock, with emissivityweighting near the nose increasing the the average shock opening angle sampled by the observer. Therefore, narrowshock opening angles the light curve exhibit a degeneracy where emissivity weighting near the shock nose can effectively hocks in MSP Binaries R for a given inclination. The different effect at larger R is due to weightingthe overall particle distribution leading to lower values of bulk motion, effectively lowering the impact of beaming, andhence the magnitude of modulation. Clearly, a self-consistent analysis is necessary for the emissivity distribution todisentangle model parameters.Finally, in Figure 9 we consider restricting θ max , X = π/
6, only considering emission from this small shock cap ratherthan the full head as in Figure 8. This scenario approximates one where emission from the shock is extraordinarilyconcentrated near the stagnation point. It is clear that such a restriction generally results in too-wide peak separationfor any inclination, ς distribution, or R values, principally due to the local flatness of this region. There is little changein light curve shape with different values of R . The emitting region is geometrically similar to large values of R (cid:38) . θ max , X = π/
2, highlighting a degeneracy in parameters for “flat” shocks. However, this degenerate regime is ruledout by the large peak separation. Therefore, the emitting region in J1023+0038 and similar IC-centered systems mustinclude synchrotron contributions from larger θ max , X values where the shock significantly curves or bows, consistentwith large radio eclipse fractions enshrouding the MSP.4.4. Discussion: Caveats and Future Refinements
The Thin-Shell Approximation
An objection to the formalism we have presented may be the use of the thin-shell approximation, especially whenknown X-ray bow shocks in PWNe do not exhibit geometrically-thin morphologies traced by radiatively coolingelectrons. Simulations, as well as observed pulsar bow shocks (e.g., Brownsberger & Romani 2014) typically show athickness of order δR /R ∼
30% near the head of the shock with increasing thickness far from the head. The thicknessis highly dependent on the physical conditions such as the in situ MHD σ (cid:46) σ (cid:28)
1, i.e., where pressure terms in the Euler equation are omittedcompared with velocity terms. For low bulk velocities near the shock nose, this approximation is clearly overstepped.However, this low-speed component near the nose is also not Doppler-boosted, barely impacting the X-ray modulationand therefore constitutes a DC or background-level offset to the X-ray modulation.Since we are principally concerned with DP light curves, only the radiative portion of the electron population,whose photons are then Doppler-boosted, are consequential in flux ratio Eq. (27). The spatial distribution of thisleptonic population is poorly understood but may be geometrically thinner than δR /R ∼ ∼ . δR /R (cid:46) . β Γ) max (cf. Figure 6). Therefore, presently, geometric thickness ofthe radiative population may be viewed as a dispensable supplemental complexity. Certainly, the formalism above isamenable to an ensemble of thin-shells with the cost of additional free parameters – a multi-zone model may embracean arbitrary assemblage of multiple single-zone thin-shells. With a single shell, quantities such as the bulk Lorentzfactor and shock density profile Σ e are to be interpreted as ensemble averages in this “one-zone” model. A naturalextension is a two-zone model, for two different electron populations separated by the contact discontinuity where oneis baryon-loaded and attains bulk velocities much lower than the relativistic adiabatic sound speed (cid:46) c/ √
3. Such anexploration is deferred for future work.4.4.2.
Diagnosing Spatially-Dependent Relativistic Particle Acceleration From Light Curves
Any spatial variation of the index p → p ( θ, φ ) induces energy dependence of the light curves in Eq. (27). Asdiscussed previously, the upstream relativistic shock geometry is quasi-perpendicular at the stagnation point relaxing4 Wadiasingh et al. to quasi-parallel at higher bulk velocity locales. Accordingly, the field compression, reconnection, the particle density,and the acceleration spectral index p probably vary from locale to locale, deviating from the canonical value of p ≈ p = 2 .
23 for ultrarelativistic parallel shocks (Kirk et al. 2000), steepening in locales whereacceleration is inefficient. Such index variation will be convolved with spatially-dependent contributions from particlecooling. The distribution ς serves as an accessible diagnostic of the spatial dependence underlying the nonthermallepton population, shock acceleration, and post-shock magnetic field. The winding of the pulsar wind in a Parker spiralmust yield different field obliquities at different angles θ along the shock surface. These obliquities will depend onthe tilt of the pulsar’s spin axis to the orbital plane, and importantly, on the plasma dynamics of its magnetosphere:the morphology of the field lines emanating from the pulsar surface is influenced by the MHD conductivity of thewind (Contopoulos et al. 2014; Kalapotharakos et al. 2014). The field obliquity angle along the shock interface canhave a critical impact upon the acceleration process. First, if the field is dynamically important, it can alter thecurvature of the shock surface as well as the level of field compression and the field orientation downstream. Next, theindex p is sensitive to the choice of the field obliquity in mildly-relativistic shocks, and the power-law index resultantfrom diffusive (Fermi-like) and shock drift acceleration is generally less than p = 2 as long as the shock interface issubluminal (Summerlin & Baring 2012). Such circumstances will likely occur at significant angles θ away from thenose or stagnation point of the shock hyperboloid. Thus, unless the microstructure of the field at the nose is highlyturbulent, the nose locale should provide a fairly steep particle distribution that flattens as θ increases, with the indexeventually declining at the lateral extremities of the shock because its MHD compression ratio will decline. Kineticturbulence will also impart a spatial dependence on the accelerated spectrum, a function of the spatially-dependentturbulent cell lengthscale (e.g., Zhdankin et al. 2016). There are clearly a number of parameters involved in describingsuch complexity.Asymmetries and energy dependence of structure in light curves also serve as a probe of the underlying particledistribution. Coriolis effects near the stagnation point break axisymmetry but are only relevant if the slow baryoniccomponent dominates the dynamics and geometry of the mildly relativistic shocked pulsar component. Observe thatinclination of the pulsar spin axis to the normal to the orbital plane may also generate an asymmetric bow shockstructure, leading to asymmetric modulation of the observed flux. Analogous to the Vela pulsar in γ -rays (Abdo et al.2010), energy dependence in ratio of the main peaks and widths may be evident. Any energy dependence of lightcurve asymmetry therefore is a signature of differential index variation rather than simple geometric asymmetry in theshock geometry.To constrain and probe such influences in the future demands high fidelity phase-resolved spectroscopy of light curves.There exists an integral equation for the flux for each observed energy bin constructed from Eqs. (19)–(23) that isfairly involved. We defer such an exploration for when such data are available, and note that such an observationalprogram on multiple systems would potentially be a remarkable astrophysical probe of particle acceleration in obliquerelativistic shocks. 4.4.3. Phase-Dependent Synchrotron Spectral Breaks and Cut-offs
A secondary source of spectral index p variation and energy dependence far beyond the classical X-ray band mayarise from transport phenomena for the steady-state particle distribution ¯ n e . Such influences could lead to a differentorigin of spectral curvature, similar in origin to the observed spectral cooling breaks in GRBs and AGNs but in anorbital-phase-dependent manner, since the selection of Doppler factors that modify the spectral elements is dependenton the observer’s viewing perspective. The underlying particle distribution and transport are constrained by wheresuch spectral breaks occur and whether they are dependent on orbital phase. Synchrotron cooling is the dominantenergy loss mechanism at large Lorentz factors γ e (cid:38) , and if the shock acceleration is efficient at the gyroscale likethe nebular shock in the Crab (de Jager et al. 1996), the radiation-reaction limited synchrotron exponential cut-off isindependent of magnetic field at ∼
160 MeV in the comoving frame; this is the rough upper limit for energy unlessextremely fast acceleration processes are occurring. Therefore, for the binary systems considered here, the convolutionof Doppler shifts from different shock locales will instill only a modest smearing in this spectral turnover in the gamma-ray waveband. Moreover, this will not impact our light curve determinations, which are germane to lower frequencies.At the other end of our broadband spectral window, simple estimates of the synchrotron self-absorption frequency(Rybicki & Lightman 1979) indicate that it is well into the radio band. This is quickly discerned for typical numberdensities (cid:104) n e (cid:105) ∼ − cm − and length scales R ∼ − cm. hocks in MSP Binaries τ conv ∼ R a/ ( βc ) in comparison to theradiative timescale τ rad in the comoving frame, τ rad = 3 m e c σ T γ e u ≈ (cid:18) γ e (cid:19) (cid:18) − u (cid:19) s , (29)where the energy density u = u rad + u B ∼ − with u B = ¯ B s / (8 π ) and u rad ≈ (cid:104) Γ (cid:105) η ˙ E SD / (4 πR c ) (efficiency η < τ conv > τ conv , min = R a/c . The radiative efficiency is, crudely, τ eff /τ rad = (1 + τ rad /τ conv ) − where τ − = τ − + τ − is the effective timescale (Tavani & Arons 1997). The radiative efficiency τ eff /τ rad (cid:28) . The inefficient cooling in the wings of the shockmay seem detrimental to the model, but observe the Doppler-boosting in flux in Eq. (21) grows faster than the declinein comoving radiative efficiency for orbital phases and shock locales where observer beaming is significant, assuming n e ∼ Q e τ eff for a spatially-independent injection rate Q e . That is, β ( θ ) in τ eff ≈ τ conv is slower growing than δ p − / D when ˆ n v · ˆ u (cid:48) ≈ /β . Alternatively, one may view it as a constraint for the unknown spatial acceleration/injection rate Q e ( γ e ), which cannot wane along the head of the shock faster than compensated by the Doppler-boosting influences.Note that the surface area element corresponding to the shock nose is also much smaller than those associated withthe wings. The fast-cooled locale is also not significantly Doppler-boosted, therefore not impactful on the background-normalized light curve morphology, but does influence cumulative spectral determinations. Accordingly, at the wingsof the shock, τ rad = τ conv > τ conv , min defines a break Lorentz factor below which radiative cooling is inefficient, whichdominates the Doppler-boosted light curves, γ e, break (cid:46) m e c R a σ T u ≈ × (cid:18) × cm R a (cid:19) (cid:18) − u (cid:19) . (30)This Lorentz factor is somewhat lower than the multi-TeV-scale γ max ∼ values expected from a purely gyroscaleacceleration radiation-reaction limited scenario in a ∼ −
100 G field. The concomitant characteristic comoving breakenergy ¯ (cid:15) c, break = 3 ¯ B s / (2 B cr ) γ e, break may be well beyond the classical soft X-ray band, (cid:15) f, break = ¯ (cid:15) f, break δ D (cid:46) m e c ) R a σ T u (cid:18) ¯ B s B cr (cid:19) δ D (31) ≈ δ D
10 MeV m e c (cid:18)
10 G¯ B s (cid:19) (cid:18) × cm R a (cid:19) , u B (cid:29) u rad . Although the characteristic comoving break energy is independent of observer perspective, the Doppler factor depen-dence of Eq. (31) introduces an orbital and geometrical influence to the break energy at any particular emission pointalong the shock. Thus, the Doppler factor must be spatially convolved with the shock emission region, introducingfactors of unity variation in the cooling break and smearing any sharp spectral transition.Observe that in the context of IC-centered systems the break energy Eq. (31) scales as r s using Eq. (1), althoughwith a σ dependence for the post-shock magnetic field B s that depends on the kinetic-scale physics in the shock. ForJ1023+0038, no phase-averaged spectral cut-off is seen by NuSTAR up to 50 −
80 keV. Assuming B s ∼ √ σB w with σ ≈ − results in r s (cid:38) cm for the shock radius constraint in J1023+0038, and requires σ < r s < a ≈ cm. During transient flaring optical states of the companion, the temporary intensification of a dominant Comptoncooling for Lorentz factors γ e (cid:46) will reduce the characteristic synchrotron break energy as u − ∼ T − . Suchspectroscopic features, phase-dependent, transient, or correlated may be observable for some bright MSP binaries withNuSTAR and a future next-generation hard X-ray to soft γ -ray Compton telescope (e.g., ComPair: Moiseev et al.2015; e-ASTROGAM: De Angelis et al. 2016). CONCLUSIONIn this paper we have focused on the radio and X-ray phenomenology of MSP binaries constructing geometricmodels for X-ray double-peaked light curves and radio eclipses. We have shown that in the thin-shell model, one may6
Wadiasingh et al. constrain the shock parameters as a function of binary inclination using radio eclipses and X-ray light curves. Thetwo different modes for the orbital phase-centering of double-peaked X-ray light curves are interpreted as owing tothe relative shock orientation, reflecting the ratio of MSP to companion wind ram pressures. When this ratio is muchlarger than unity, orbital radio eclipse fractions of the MSP f E may be low and X-ray DP light curves are centered atsuperior conjunction of the pulsar. In contrast, when this ratio is small, the phase-centering is at inferior conjunctionand f E is large. This shock orientation along with large f E advances the scenario where the pulsar is enshrouded bythe shock somewhere past the L point. Synthetic X-ray light curves are a complex function of orbital inclination,shock geometry, bulk Lorentz factor, and emissivity distribution. In the power-law regime, the light curve morphologyis energy-dependent if and only if the power-law index of the steady-state particle distribution is spatially dependentalong the shock. This then is a powerful probe of particle acceleration in relativistic oblique shocks.We have argued that the shock in BWs and RBs have two components that are not well-mixed unless the companionmass loss rate is low, consistent with the reasoning that relatively low baryon loading is necessary for the mildlyrelativistic 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 stableshock at orbital length scales in our model. Accordingly, some mechanism is necessary to stabilize the shock fromgravitational 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-energyemission and transport processes. The quantitative agreement of the synthetic light curves with observations stronglyencourages further development of the model to include considerations of diffusive and convective transport in the shockenvirons. 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 thenear future of DP X-ray light curves may be able constrain the shock physics given orbital parameters, or constrainorbital parameters such as binary inclination (and consequently MSP mass) using a model for the shock geometry andthe particle distribution in different regions. Energy-dependent light curves by NICER or NuSTAR may also probe theunderlying relativistic shock acceleration’s spatial dependence in the future. Our study lastly motivates populationstudies of radio eclipse phenomenology in the MeerKAT/SKA era, where the number of MSPs discovered will increaseby severalfold (Keane et al. 2015).We thank the anonymous referee for aiding to improve the flow and structure of the manuscript. Z.W. acknowledgeshelpful discussions with Cees Bassa, Julia Deneva, Guillaume Dubus, Jason Hessels, Christopher Johns-Krull, PatrickKilian, Edison Liang, Alessandro Papitto, Martin Pohl, Scott Ransom, Mallory Roberts, Bronek Rudak, Ben Stappers,and Kent Wood. C.V. and Z.W. acknowledge Tunde Ayorinde’s aid in some cross-checks. C.V. & Z.W. are supportedby the South African National Research Foundation (NRF). The work of M.B. is supported by the South AfricanResearch Chairs Initiative of the Department of Science and Technology and the NRF. Any opinion, finding andconclusion or recommendation expressed in this material is that of the authors and the NRF does not accept anyliability 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 SYSTEMSA.1.
General Formalism
Although the focus of this text is on eclipses by optically-thick intrabinary shocks, the analytical formalism herecan also be applied to arbitrary azimuthally symmetric surfaces that occult a point source. The case of eclipses byquasi-static Roche lobes has been treated analytically in Chanan et al. (1976) and Kopal (1959), however the methodpresented 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 ˆ z directionwith binary angular frequency Ω b . The observer angle is defined such that ˆ z · ˆ n v = cos i with ˆ n v = cos i ˆ z + sin i ˆ x theobserver line-of-sight unit vector. To obtain the projection of the binary system into the plane of the sky perpendicularto ˆ n v , we construct a new primed rotated coordinate system such that ˆ x (cid:48) is parallel to ˆ n v at arbitrary orbital phase.If the cartesian basis ˆ r → (ˆ x , ˆ y , ˆ z ) (cid:124) defines the vector space at phase Ω b t = 0 and inclination i = π/
2, then the primed hocks in MSP Binaries i = 0 and i = π/ b t = 0 prescribing the superior conjunction where the companion is between the MSP andobserver in Eq. (A3) of the orbital equations of motion. The rotations are given by ˆ r (cid:48) = Λ i Λ Ω b t ˆ r , (A1)where Λ Ω b t = cos (Ω b t ) − sin (Ω b t ) 0sin (Ω b t ) cos (Ω b t ) 00 0 1 Λ i = sin i i − cos i i (A2)define the primed coordinate basis about the barycenter such that ˆ x (cid:48) = ˆ n v , with the span of ˆ y (cid:48) and ˆ z (cid:48) characterizingthe plane of the sky.The stars are treated as following unperturbed Keplerian trajectories, although a post-Newtonian generalization isstraightforward, with phase zero constructed to rest on the ˆ x -axis with the following equations of motion,MSP (primary): r NS ( t ) = − r NS ( t ) [cos (Ω b t ) ˆ x + sin (Ω b t ) ˆ y ] (A3)Companion (secondary): r c ( t ) = r c ( t ) [cos (Ω b t ) ˆ x + sin (Ω b t ) ˆ y ] . The scalars r i ( t ) for non-zero eccentricity depend on the orbital phase, r i ( t ) = r i (1 − (cid:15) )1 + (cid:15) cos (Ω b t ) , r i = r c , r NS (A4)where r c + r NS ≡ a ) are the respective semi-major axes for each star, with r c = q/ ( q + 1) and r NS = 1 / ( q + 1) where q = M MSP /M c is the mass ratio. For nonzero eccentricity, which is beyond the scope of thispaper, 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 secondaryare given by r (cid:48) NS ( t ) = r NS ( t ) (cid:104) − sin i cos (Ω b t ) ˆ n v − sin (Ω b t ) ˆ y (cid:48) + cos i cos (Ω b t ) ˆ z (cid:48) (cid:105) (A5) r (cid:48) c ( t ) = r c ( t ) (cid:104) sin i cos (Ω b t ) ˆ n v + sin (Ω b t ) ˆ y (cid:48) − cos i cos (Ω b t ) ˆ z (cid:48) (cid:105) with a parallel/orthographic projection on the Span { ˆ y (cid:48) , ˆ z (cid:48) } plane comprising the observer view.For orbital angular speeds a Ω b (cid:28) c , the coordinates of a vector (cid:96) defining a ray leaving the primary towards theobserver traversing a distance l (in units of a ) can be expressed as (cid:96) (cid:48) = r (cid:48) NS ( t ) + l ˆ n v , (A6)with the distance between an interaction point at the ray and an arbitary location in the system r (cid:48) m (e.g., thesecondary) given by (cid:12)(cid:12) (cid:96) (cid:48) − r (cid:48) m (cid:12)(cid:12) . The optical depth for a specified absorption coefficient α is computable in the usualway through the observer line-of-sight integral at each orbital phase, τ = (cid:82) αdl . Such a procedure for pulsar eclipsesby 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 iscomputationally inefficient for arbitrary volumes when the optical depth is large for the region of interest, and whenthe transition region of low-to-high optical depth is sharp. In the case of optically thick surfaces and volumes, weemploy a geometric occultation formalism described below.For a surface with azimuthal symmetry with radial function R ( θ ), every locus of points at fixed θ is a circle. Whenprojected onto the Span { ˆ y (cid:48) , ˆ z (cid:48) } plane of the sky, these circles are transformed to ellipses by the elementary rotationsof orbital phase and inclination angle. The definition of eclipses by a bow shock is thus transformed to the problemof determining if the projected location of the MSP is interior to any projected bow shock ellipse, for all prescribed8 Wadiasingh et al. θ . We define θ max , R as the maximum polar angle where the shock surface transitions from optically-thick to thin ata given frequency. The maximum length along the axis of symmetry of the shock tail past the companion position L z ( θ max , R ) = − R ( θ max , R ) cos θ max , R should be smaller than the typical orbital length scale a for flow velocities thatare similar in scale to the orbital speed. These relatively slow flow velocities are interpreted to originate from thenonrelativistic baryonic component of the companion wind, rather than the leptonic pulsar wind, with the two windsseparated by a contact discontinuity.The locus of points for a particular projected ellipse, given the conventions developed above, can be found to be E ( t, θ, φ ) = ( α y + κ y sin φ ) ˆ y (cid:48) + ( α z + β z cos φ + κ z sin φ ) ˆ z (cid:48) , (A7)with α y = S ( t, θ ) sin (Ω b t ) α z = − S ( t, θ ) cos (Ω b t ) cos iβ z = R ( θ ) sin θ sin i (A8) κ y = R ( θ ) sin θ cos (Ω b t ) κ z = R ( θ ) sin θ sin (Ω b t ) , where φ ∈ (0 , π ) is the azimuthal angle parameterizing the azimuthally-symmetric surface. The function S ( t, θ ) defineswhether the shock surface surrounds the primary or secondary with the orbital phase convention defined above: forsurrounding (bowing around) the secondary it takes the form S c = r c ( t ) − R ( θ ) cos θ but S NS = − r NS ( t ) + R ( θ ) cos θ for where the shock enshrouds the pulsar.With these definitions, it is evident that the center of the ellipse is then { α y , α z } . In general, the ellipses’ axes ofsymmetry are rotated by an angle with respect to the ˆ y (cid:48) or ˆ z (cid:48) axes. The semi-major E a and semi-minor E b axeslengths are easily shown to be E a = | R ( θ ) sin θ | (A9) E b = | R ( θ ) sin θ sin i cos(Ω b t ) | , for a given inclination and orbital phase and, for a given θ . The projected eclipsed area is ∝ R ( θ ) sin θ sin i cos(Ω b t ).Using elementary methods, the angle between an ellipse’s semi-major axis and the ˆ y (cid:48) direction is found to becos ϑ E = 2 | cos (Ω b t ) | cos i (cid:113) i ) − b t ) sin i , (A10)for phases − π/ ≤ Ω b t ≤ π/ G x ( t ) = (cid:104) r (cid:48) NS · ˆ y (cid:48) − α y (cid:105) cos ϑ E + (cid:104) r (cid:48) NS · ˆ z (cid:48) − α z (cid:105) sin ϑ E (A11) G y ( t ) = − (cid:104) r (cid:48) NS · ˆ y (cid:48) − α y (cid:105) sin ϑ E + (cid:104) r (cid:48) NS · ˆ z (cid:48) − α z (cid:105) cos ϑ E , which can be used to test if the projected position of the MSP is eclipsed. We define a unit Heaviside test functionΘ(Γ) with argument Γ = (cid:18) G x E a (cid:19) + (cid:18) G y E b (cid:19) − , (A12)is less than zero for eclipses. Hence the position G x , y and parameter Γ can be associated with a physical attenuation(and τ ) at any given θ , i.e., e − τ ∼ f ( G x , G y ) for some function f ; we choose f to be a unit step function in ouroptically thick formalism. Then, Θ(Γ) is unity when the pulsar is not eclipsed. The eclipsing by the entire ensembleof ellipses j max that encompass the employed surface is a product of Θ for every θ up to θ max , R , suitably discretizedto sample the complete surface, Θ tot = j max (cid:89) j Θ[Γ( θ j )] . (A13) hocks in MSP Binaries θ max , R . This yields a lowerlimit because a general surface may have pockets of plasma that may occlude the emission source in a nontrivialmanner. For cases where the bow shock is not swept back and the curvature is positive everywhere, the precedingconstruction is unnecessary, since the shrouding fraction is analytically expressible by Eq. (13).It is numerically expedient to adopt a continuous approximation of the Heaviside function, Θ[Γ( θ j )] ≈ / (1 + e −T Γ( θ j ) ), where T is an associated pseudo optical depth that parametrizes the geometric transition thickness nearthe occluding surface R ( θ ) and approximates a sharp transition for T (cid:29) tot ≈ exp − j max (cid:88) j log (cid:0) e −T Γ j (cid:1) T (cid:29) . (A14)Because of the symmetry in the problem, the full orbit and other branches for the solution of cos ϑ E need not beconsidered, and the calculation can be restricted to one quarter of the total orbit. The total eclipse fraction f E by thesurface during a total duration 2 π of an orbit is then given by the integral, f E = 1 π (cid:90) π/ d (Ω b t ) [1 − Θ tot (Ω b t )] (A15)for a given set of Keplerian orbital parameters and radial function R ( θ ).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 ateach instantaneous orbital phase then it is straightforward to include the sweep-back and Coriolis forces due to orbitalmotion by retarding the surface at each θ by a specified time delay. We neglect Coriolis effects perpendicular to theoutward direction, equivalent to the assumption that v ⊥ (cid:28) v inj where v ⊥ is the bulk flow velocity perpendicular tothe line connecting the two stars; this approximation, which preserves axisymmetry of the shock, is generally valid forlocales distant from the stagnation point, the principal focus of this work associated with radio eclipses. Connectingorbital phases to the surface parameter θ then suffices to produce a good approximation of such swept-back shocktails, and the eclipse fraction formalism of the previous section can be expediently utilized. This ballistic assumptionis valid if the shock bulk velocity is highly supersonic, as in MSP binaries, especially if the winding number of theswept-back tail that contributes to the eclipses is zero. The assumption v ⊥ (cid:28) v inj requires R (cid:28) . θ (cid:38) π/ θ or longer tails. Therefore, the validity ofthis approach is restricted to the parallel-wind geometry explored in Appendix B.The model we adopt here simply injects the shock flow at a finite velocity v inj near the companion at θ = π/ v inj (cid:29) v esc ∼ cm s − which, coincidently, is also the same orderof magnitude as the orbital velocity v 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 wind (in the lab nonrotating frame) dueto radiation pressure from the pulsar wind is δt = − v inj + (cid:113) a wind aL z ( θ ) + v a wind , (A16)which reduces to aL z /v inj for small flow accelerations. For zero acceleration, the approximate tail angle with respectto the line connecting the two stars is then arctan( v orb /v inj ), good when v inj (cid:29) v orb . For large a wind or v inj , the timedelay approaches zero and the symmetric shock solution is recovered. Here L z ( θ ) = − R ( θ ) cos θ connects the timedelay to points along the shock surface. Thus, generically, increasing asymmetry increases the total eclipse fraction ata given θ max , R .The computation of eclipse fraction is straightforward with the simple replacement Ω b t → Ω b ( t (cid:48) − δt ) in Eq. (A12)at each θ > π/ L z or θ max , R for all variables associated with the shock surface. The ingress0 Wadiasingh et al.
Figure 10.
A similar suite of axisymmetric computations as that of Figure 2 but for the parallel-wind shock geometry, withdifferent but appropriate choices θ max , R in the rightmost two panels corresponding to shock tail lengths past the companion of R and 4 R . and egress of eclipses are defined to occur for phases Ω b t (cid:48) ∈ [ − π/ ,
0) and Ω b t (cid:48) ∈ (0 , π/
2] respectively. The totaleclipse fraction is then f E = f E (Ingress) + f E (Egress), and the asymmetry can be characterized by the difference∆ f E ≡ f E (Egress) − f E (Ingress) and the ratio of egress to ingress eclipse durations. B. PARALLEL WIND BOW SHOCKAn analytic alternative to the colliding isotropic winds shock geometry of Eq. (6) is that of a parallel-isotropic windinteraction, typically invoked in the context of bow shock nebulae (Wilkin 1996). The radial function, defined on θ ∈ (0 , π ), is manifestly scale-invariant R para ( θ ) R = csc θ (cid:112) − θ cot θ ) . (B17)The geometry also constitutes the R (cid:28) R ≈ η w (cid:28) θ → π .This geometry may be more relevant to SC-centered BWs such as B1957+20 where the companion pressure supportfrom an irradiated evaporative wind is anisotropic and weaker on the night side. The net effect of such anisotropieswould narrow the shock opening angle versus that of Eq. (6) for the same stagnation point R value. The two geometriesconverge only for R values well below the companion radius, R ∗ ∼ .
1. Therefore for R ≥ R ∗ , the geometry Eq (B17)is much narrower than that of Eq. (6). Moreover, the tail region θ (cid:38) π/ θ , in dramaticcontrast to the isotropic-winds scenario; the ratio of perpendicular to longitudinal shock components (with respect tothe shock symmetry axis) tends to zero as θ → π .B.1. B1957+20 Radio Eclipses
When the shock enshrouds the companion, scale-invariance of Eq. (B17) is broken for eclipses and f E is dependenton R , as for the two isotropic-winds geometry in Figures 2–4. We compute the parallel-wind geometry counterpartsof these Figures as applied to B1957+20, mirroring § θ max , R left-to-right.However, unlike its counterpart Figure 2, the two rightmost panels prescribe values of θ max , R independent of R owingto the scale-invariance of Eq. (B17). The values θ max , R ≈ { . , . } , independent of R , prescribe shocks with tails oflongitudinal lengths | R ( θ ) cos θ | = L z = { R , R } . The narrowness of the shock geometry also imparts a much slowergrowth of f E as a function of R , and implies larger values of R for the same f E than the isotropic-winds case. As forthe isotropic-winds case, there are constraints on θ max , R . For i = 85 ◦ , the implication is unchanged: θ max , R (cid:46) π/ R > R ∗ . Observe that for i = 65 ◦ and f E (cid:46) θ max , R > { π/ , π/ , . , . } corresponds to limits R (cid:46) { . , . , . , . } and is of lower sensitivity to the value of θ max , R than the isotropic-winds geometry. Theorigin of this behavior is readily apparent from inspection of the growth curves in Figure 11 – there is a plateauing of f E for large values of θ max , R , following a steep rise through the head of the shock. This result intrinsic to the geometry hocks in MSP Binaries Figure 11.
Eclipse fraction growth rates for the parallel-wind shock geometry, contrasting the isotropic-winds case in Figure 4.The case i = 55 ◦ requires R > .
35 for non-zero f E for all θ max , R and is omitted. and owing to narrowing of the shock in the tail. Low inclinations are therefore untenable with this shock geometry,with i (cid:46) ◦ requiring R values that may exceed 0 . f E (cid:38) § A.2 we may explore the asymmetry of eclipses for the long tails of theparallel-wind geometry. We assume this asymmetry arises from the relatively slow and ablated shocked companionwind with eclipsing locales being far behind the companion in the downstream tail of the shocked companion flow.The shock is assumed to be azimuthally symmetric at every instantaneous orbital phase at the companion position,an expedient approximation that is accurate when the transverse velocity component of the shock is much smallerthan the parallel component v ⊥ (cid:28) v inj , which exceeds the escape velocity v inj (cid:38) v esc ∼ cm s − . To simulate thetail sweepback, a supersonic flow of velocity v inj (cid:29) v esc is prescribed for the downstream shocked companion wind,accelerated by radiation pressure from the pulsar wind by parameter a wind parallel to the line joining the two stars,neglecting gravitational forces. The flow must be supersonic in this ballistic model, so that the internal hydrodynamicsof the flow can be neglected, e.g., Rasio et al. (1989), and we go beyond previous analyses by considering generalinclination angles. Using the analytics developed in § A.2, the ingress and egress eclipse fractions can be computedfor a given R , inclination i , tail length L z measured from the companion position, v inj , and a wind . In Figure 12, wedepict two different cases of constant velocity and constant acceleration, and it is evident that acceleration of the flowcan have dramatic consequences on the geometry and eclipse asymmetry depending on the parameters L z , a wind , and v inj , which unfortunately, are poorly constrained observationally due to degeneracies in these parameters for a givenasymmetric eclipse.From observations of B1957+20 in Ryba & Taylor (1991), the total eclipse fraction is constrained to be conservatively0 . (cid:46) f E (cid:46) .
15, while the egress minus ingress eclipse difference about SC, 0 (cid:46) ∆ f E (cid:46) .
05, with largest asymmetryfor the lowest observing frequencies. In Figures 13 and 14 we explore the parameter space for B1957+20 computing∆ f E as a function of v inj versus L z = | R ( θ ) cos θ | or a wind , respectively, excluding portions of the parameter spacewhere the total eclipse fraction or asymmetry is disallowed by observations. In these figures, R is modified uponchanging i to keep the total eclipse fraction similar. For fixed L z , the asymmetry decreases with increasing v inj and theflattening in the asymmetry at larger L z is due to the long shock tail. If the spatial steady-state density monotonicallydecreases with L z downstream, then larger distances naturally produce larger eclipse asymmetry sampled at lowerobserving frequencies. Horizontal cuts in Figure 13 trace the frequency dependence of the eclipse asymmetry in2 Wadiasingh et al.
Figure 12.
The effect of constant velocity and acceleration on the parallel-wind swept-back shock tail, out to unity semi-major distance a , in the PSR B1957+20 system due to orbital motion and Coriolis effects. The system angular momentumvector direction is pointing out of the plane in the top two panels, i.e. the system is rotating counterclockwise. The depictedorbital phase is at eclipse egress. For all panels, the dark green coloring mimics the loci of points corresponding to the contactdiscontinuity injected supersonically at velocity v inj = 10 cm s − at the companion position, while the brown colored linesillustrate a symmetric unswept high-speed shock component. The left panels highlight the case of zero wind acceleration, whilethe right panels accelerate the shock at a constant rate a wind = 10 . cm s − (outward, parallel to the line joining the two stars).The constant acceleration, contrasted to zero acceleration, decreases the total eclipse fraction from 11 .
33% to 8 . f E = 4 .
75% to 1 .
9% and egress-ingress ratio from 2 .
45 to 1 .
58. The purple points denote the geometriccenter of the projected ellipses from Eq. (A8). this simplified model of eclipse mapping. For a modest tail length of L z ∼ v inj (cid:38) cm s − although smaller values approaching the escape speed ∼ cm s − are allowed if theacceleration a wind (cid:38) cm s − . The allowed parameter space is generally larger for lower i . There are obviously hocks in MSP Binaries Figure 13.
The asymmetry of eclipses given by the proxy ∆ f E for the B1957+20 system due to orbital sweep-back of aparallel-wind shock, as a function of injected velocity, with a wind = 0, at the companion and shock length L z = | R ( θ ) cos θ | downstream from the companion to a maximum of 4–5 R . The white region corresponds to the portion of the parameter spacenot allowed by the radio results of Ryba & Taylor (1991); here either the asymmetry of eclipses or the total eclipse fractionare too large even at the lowest frequencies. The orbital speed of the companion is shown as the blue line and varies withinclination due to observed pulsar mass function for a fixed mass ratio q = 69 .
2. Larger speeds or shorter tails reduce totaleclipse asymmetry. Horizontal cuts at a fixed v inj in the figure may be interpreted as the frequency dependence of the eclipseasymmetry, with lower frequencies sampling larger L z . Figure 14.
The asymmetry of eclipses given by the proxy ∆ f E for the B1957+20 system due to orbital sweep-back of a parallel-wind shock, as a function of injected velocity and constant wind acceleration a wind for a fixed shock length L z = 1 downstreamfrom the companion. Horizontal cuts depict the eclipse asymmetry ∆ f E varying the a wind at a fixed injection velocity, whilevertical cuts correspond to ∆ f E as a function of v inj . The white region corresponds to the portion of the parameter space notallowed by the radio results of Ryba & Taylor (1991); here either the asymmetry of eclipses or the total eclipse fraction aretoo large even at the lowest frequencies. Larger speeds and accelerations reduce eclipse asymmetry. The orbital speed of thecompanion is shown as the blue line and varies with inclination due to observed pulsar mass function for a fixed mass ratio q = 69 . Wadiasingh et al. % β max = i = o , R = ϛ ( θ ) i = o , R = ϛ ( θ ) i = o , R = ∂ϛ ( θ ) ∂θ = i = o , R = ∂ϛ ( θ ) ∂θ = i = o , R = R * , ϛ ( θ ) i = o , R = R * , ∂ϛ ( θ ) ∂θ = - π - π π π Angular Orbital Phase R e l a t i v e L u m i n o s i t y % β max = i = o , R = ϛ ( θ ) i = o , R = ϛ ( θ ) i = o , R = ∂ϛ ( θ ) ∂θ = i = o , R = ∂ϛ ( θ ) ∂θ = i = o , R = R * , ϛ ( θ ) i = o , R = R * , ∂ϛ ( θ ) ∂θ = - π - π π π Angular Orbital Phase R e l a t i v e L u m i n o s i t y Figure 15.
Parallel-wind shock geometry synchrotron light curves contrasting those in Figure 7. Note the different value of R for i = 65 ◦ required for consistency with radio eclipses. several confounding variables and large degeneracies in the parameter space for such non-axisymmetric eclipses at afixed observational frequency that may produce different combinations of f E and ∆ f E . This motivates radio populationstudies of eclipsing MSP binaries and eclipse phenomenology to unearth commonalities and constrain the parameterspace of the two-wind interaction. B.2. B1957+20 Synthetic X-ray Light Curves
Mirroring Figure 7, we compute light curves using the parallel-wind geometry for the head of the shock θ max , X = π/ R = { . , . } consistent with eclipses for the head of the shock at i = { ◦ , ◦ } ,respectively, as well as the limiting case R = R ∗ . Unlike for the isotropic-winds geometry, the value of R in the scale-invariant parallel-wind geometry does not regulate peak separation or width, only the influence of shadowing governedby the scale R ∗ ≈ .